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Exact  conditions  for  the  existence  of  limit  cycles  due 
to  friction  in  a  servomechanism  have  been  determined  for  the 
general  (n-dimensional)  case.   These  conditions  are  in  the 
form  of  nonlinear  algebraic  equations,  for  which  a  solution 
consistent  with  certain  assumptions  must  be  found. 

A  two-input  model  for  friction  is  required  for  accurate 
results,  including  both  dry  friction  and  static  friction 
components.   Therefore,  standard  nonlinear  analysis  methods 
such  as  absolute  stability  and  describing  functions  have 
limited  applicability.   A  piecewise  linear  method  is  used, 
which  constructs  the  limit  cycle  trajectory  with  pieces  from 
each  region  of  the  state  space. 

The  system  model  is  first  converted  to  a  control 
canonical  form,  which  simplifies  the  analysis.   Assumptions 
made  in  the  formation  of  the  model  are  specified  exactly; 
the  model  represents  a  rotational  servomechanism,  but  is 
shown  to  be  applicable  to  other  systems. 

vi 


Next,  the  piecewise  linear  method  is  applied  to  an 
assumed  limit  cycle.   The  resulting  algebraic  equations  can 
be  solved  for  the  exact  limit  cycle  period,  amplitude,  and 
trajectory,  if  any  exist.   The  limit  cycles  found  are  of 
simple  type;  that  is,  complex  orbits  such  as  those  which 
traverse  a  region  more  than  once  would  not  be  found. 
Necessity  and  sufficiency  of  the  existence  conditions  is 
proved. 

An  analysis  of  the  local  orbital  stability  of  limit 
cycles  found  by  this  method  is  also  presented.   Eigenvalues 
of  a  given  matrix  determine  the  asymptotic  stability  of  the 
periodic  orbit. 

The  property  of  symmetry  of  limit  cycles  for  odd 
systems  in  general  and  friction  in  particular  is  examined. 
Special  cases  are  presented  in  which  symmetry  holds,  and 
counter-examples  for  the  general  case  of  odd  systems  are 
given.   The  question  of  symmetry  for  friction  limit  cycles 
remains  open. 

Examples  are  included  to  illustrate  the  method.   One 
example  is  a  3D  system  that  exhibits  two  limit  cycles  in  its 
phase  portrait,  each  a  different  type. 

Listings  of  useful  analysis  computer  programs  are 
given,  in  addition  to  suggestions  for  further  research  in 
this  area. 
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CHAPTER  I 
INTRODUCTION 

Problem  Statement 

The  objective  of  this  research  is  to  define  exact 
conditions  on  the  existence  and  behavior  of  limit  cycles  due 
to  friction.   The  analysis  method  to  be  used  is  a  piecewise 
linear  approach  which,  unlike  approximate  methods  such  as 
describing  functions  or  absolute  stability,  has  the 
potential  for  providing  exact  (i.e.,  necessary  and 
sufficient)  conditions  for  existence. 

This  problem  is  of  interest  in  the  design  of 
servomechanisms,  which  must  deal  with  the  nonlinearities  in 
typical  electromechanical  devices,  such  as  friction. 
Experience  has  shown  that  the  effect  of  friction  is 
generally  benign,  since  it  tends  to  provide  increased 
damping  to  the  system.   However,  this  cannot  be  guaranteed 
for  all  systems.   Furthermore,  it  is  a  standard  practice  to 
precompensate  for  this  effect  by  deliberately  tuning  the 
system  to  be  underdamped  without  friction,  so  that  it  is  not 
excessively  damped  in  actual  operation  with  friction.   It  is 
postulated  that  this  technique,  when  pushed  to  extremes, 
results  in  large-signal  instability,  since  the  damping 
effect  of  friction  varies  inversely  with  amplitude.   It  is, 
therefore,  a  significant  research  problem  to  determine  the 
point  v/here  limit  cycles  can  exist. 
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An  approach  has  been  developed  during  this  research 
that  uses  the  piecewise  linear  method,  that  provides  an 
exact  condition  for  the  existence  of  a  limit  cycle  in  terms 
of  a  nonlinear  algebraic  system  of  equations.   When  a 
solution  to  the  system  is  found  (and  thus  a  limit  cycle 
exists) ,  the  method  automatically  gives  exact  frequency, 
amplitude,  and  state  trajectories  for  the  limit  cycle.   On 
the  other  hand,  exact  conditions  under  which  the  nonlinear 
system  of  equations  is  solvable  are  more  difficult  because 
of  the  nonlinear  nature  of  the  system. 

The  specific  problem  under  study  is  shown  in  block 
diagram  form  in  Figure  I-l,   The  plant  in  the  standard 
system  is  represented  by  a  Laplace-transform  model  of  a 
motor/load  combination,  which  includes  the  effect  of  load 
inertia  and  damping.   The  model  does  not  include  compliance 
in  the  load,  although,  as  described  in  Chapter  III,  this  may 
be  included.   The  inertia  (J)  and  damping  (B)  represent  the 
lumped  quantities  for  the  load  and  motor,  with  any  gear 
ratio  taken  into  account.   The  servo  controller  is 
represented  as  a  linear  system  with  motor  drive  torque  as 
the  output  quantity.   Note  that  despite  the  specific 
description  of  the  problem  model,  it  will  be  shown  in 
Chapter  III  (Development  of  Standard  Model)  that  the  model 
is  actually  quite  general. 

The  model  used  for  nonlinear  friction  includes  both 
static  and  sliding  coulomb  friction  components  (see  Figure 
1-2) .   Therefore,  it  is  a  two-input  nonlinearity  and  is 


INPUT  =  0 


CONTROLLER 
G(s) 


FIGURE  I-l.   BLOCK  DIAGRAM  OF  SYSTEM  UNDER  STUDY 


SLIDING  MODE 

(FOR  NONZERO  VELOCITY,  Y^;^0) 

L,=    L^sgn[Y^] 


L    (TORQUE) 


(VELOCITY) 


STATIC  OR  SUCTION  MODE 
(FOR  ZERO  VELOCITY,  Y  =0) 

Lf=     Lssat[L,^/Ls] 


L    (TORQUE) 


(INPUT 
TORQUE) 


FIGURE    1-2.       TWO-INPUT    FRICTION    MODEL 


4 
thus  difficult  to  handle  by  standard  techniques  (such  as 
describing  functions  or  absolute  stability  criteria) .   The 
extra  complexity  caused  by  the  use  of  the  two-input  model  is 
necessary;  examples  are  known  where  modelling  sliding 
friction  only  indicates  a  limit  cycle  exists,  when,  in  fact, 
the  "sticky"  effect  of  the  static  friction  prevents  it.   On 
the  other  hand,  the  sliding  friction  can  provide  much  of  the 
damping  of  many  servomechanisms  which  have  little  or  no 
viscous  friction  damping  in  the  load.   These  com.ments 
indicate  that  a  model  including  both  effects  is  necessary. 
Note  that  the  friction  model  used  is  odd,  and  therefore  the 
system  derivative  vector  will  be  an  odd  function  of  the 
state. 

The  Piecewise  Linear  Method 
The  piecewise  linear  method  is  an  approach  where 
systems  with  piecewise  linear  nonlinearities  can  be  analyzed 
using  the  tools  of  linear  system  theory.   Since  the 
nonlinearity  is  linear  over  portions  of  its  domain,  the 
state  space  can  be  divided  into  regions  in  which  the  system 
behaves  in  a  purely  linear  fashion.   The  linear  system 
theory  then  applies,  although  there  will  be  a  different 
linear  system  in  each  region.   Solution  trajectories  must 
then  be  pasted  together  from  different  regions  through  the 
use  of  boundary  conditions  (called  commutation  equations  in 
relay  servo  or  on-off  system  analysis  such  as  in  Weischedel, 
1973)  . 


The  nonlinear  friction  to  be  analyzed  in  this  study  is 
piecewise  linear,  even  though  it  requires  two  inputs.   Five 
regions  of  the  state  space  are  required  (Figure  1-3  shows 
the  case  for  2-dimensional  state  space) : 

(I)  velocity  <  0, 

(II)  velocity  >  0, 

(III)  velocity  =  0  and  magnitude  of  accelerating  input 
torque  <  Lg , 

(IV)  velocity  =  0  and  torque  >  Lg ,  and 

(V)  velocity  =  0  and  torque  <  -Lg, 

where  Lg  is  the  maximum  static  friction  torque  (breakaway 
torque) .   The  first  two  regions  are  always  open  half-spaces, 
while  the  last  three  divide  the  (n-1) -dimensional  sub-space 
where  velocity  equals  zero  into  three  regions.   For  example, 
for  the  3D  case  (n=3) ,  the  last  three  regions  are  portions 
of  the  plane  y3  =  velocity  =0.   As  can  be  seen  by 
examination  of  the  friction  model,  the  nonlinearity  either 
has  a  constant  output  throughout  each  region  (I,  II,  IV,  and 
V) ,  or  its  output  is  proportional  to  a  linear  combination  of 
states  (region  III) .   Even  a  simple  limit  cycle  must 
traverse  at  least  three  of  these  regions,  so  the  linear 
systems  of  each  must  be  considered  in  this  limit  cycle 
analysis. 

There  are  several  advantages  in  the  use  of  the 
piecewise  linear  technique.   Once  a  limit  cycle  is  found 
(i.e.,  solutions  from  four  regions  that  can  be  pasted 
together  to  form  a  closed  trajectory) ,  the  exact  trajectory 


LINEOF_ 
ZERO 

INPUTTORQUE  "n 
Y2  =  -(K/B)  Y1 
(FOR  2D  CASE) 


REGION  IV: 
VELOCITY  =  0 
ANDTORQUE>LS 


VELOCrPi' 
Y2  OR  DY1/DT 


REGION  II: 
VELOCITY  >  0 


SAMPLE  PERIODIC 
ORBIT 


REGION  III: 

UNE  SEGMENT 
OF  EQUIL  POINTS 

LS  /  K   <  Y1  <    LS  /  K 


POSITION 


Y1 


REGION  V: 

VELOCITY  =  0 

AND  TORQUE  <-LS 


REGION  I: 
VELOCITY  <0 


FIGURE    1-3.       PIECEWISE-LINEAR    REGIONS    FOR    2D    CASE 


7 
is  known.   Therefore,  the  limit  cycle  amplitude,  period,  and 
characteristics  are  known  exactly,  rather  than  approximately 
as  in  describing  functions.   Furthermore,  the  method  leads 
to  a  more  exact  analysis,  since  the  commutation  conditions 
require  the  solution  trajectories  to  match  exactly  at  the 
region  boundaries.   Finally,  the  technique  has  a  geometric 
approach  that  is  intuitively  appealing  and  lends  itself  to 
graphical  solution. 

The  method  is  applicable  to  the  majority  of  servo 
nonlinearities.   Backlash,  friction,  and  most  other 
nonlinearities  found  in  practice  are  piecewise  linear. 

The  disadvantage  of  the  method  is  in  the  difficulty  of 
solving  the  commutation  conditions  for  a  solution--f inding 
the  trajectories  that  will  "paste"  together  to  form  a  cycle. 
One  way  to  surmount  this  problem  applies  the  method  as  a 
secondary  step  in  the  analysis  in  order  to  refine  an 
approximate  solution  found  by  other  methods.   For  example,  a 
limit  cycle  period  and  amplitude  are  determined 
approximately  by  describing  functions;  then  the  piecewise 
linear  method  is  used  (with  the  approximate  solution  as  an 
initial  solution  trajectory)  in  an  iterative  way  to  converge 
to  the  true  solution.   A  composite  computer  analysis  program 
could  integrate  these  functions  for  the  convenience  of  the 
servo  designer.   Therefore,  the  difficulty  of  solution  does 
not  really  limit  the  applicability  of  the  method. 
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Illustrative  Examples 
Before  going  into  the  detailed  development  of  the 
method  and  other  results,  it  is  useful  to  present  simple 
example  systems  exhibiting  friction  limit  cycles.   The  limit 
cycles  were  analyzed  by  the  methods  to  be  presented  later, 
so  we  can  return  to  these  examples  at  the  appropriate  time 
to  illustrate  the  method;  here  only  the  results  are  shown. 
The  second  example  system  is  particularly  interesting 
since  it  is  only  3D  (three  state  variables  describe  the 
system) ,  yet  it  exhibits  two  distinct  limit  cycles  of 
differing  types  (when  the  appropriate  parameter  values  are 
selected) .   The  fact  that  such  a  simple  system  can  have  such 
complex  behavior  indicates  the  richness  of  the  study  of 
piecewise  linear  systems. 

Example  I-l;   Two-Dimensional  (2D)  System  with  Sliding 
Friction  Limit  Cycle 

The  question  of  friction  limit  cycles  for  the  case  n=2 
can  be  completely  solved;  in  fact,  the  piecewise  linear 
method  is  the  same  as  standard  phase  plane  methods  for  this 
case.   For  n=2 ,  the  controller  can  be  a  constant  gain  only 
(since  there  are  two  states  in  the  plant) ;  the  state 
variable  model  of  the  system  is  then 
(1.1)      dyi/dt   =  y2 

dy2/dt   =   -(K/J)  yi   -   (B/J)  ys   -  L^/J 
where  Lp  is  the  friction  torque  from  the  model  in  Figure 
1-2,  and  any  feedback  gain  on  y2  is  lumped  in  with  B.   It 
can  be  shown  (see  Thaler  and  Pastel,  1962,  pp.  96-104)  that 
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for  B>0  and  for  K>0,  the  system  is  globally  asymptotically 
stable  (i.e.,  no  limit  cycles).   Thus  a  stable  system 
remains  stable  (in  fact  it  is  more  stable)  with  friction  (2D 
case  only) . 

The  case  B=0  is  also  stable;  the  linear  system  has 
imaginary  eigenvalues,  the  friction  gives  enough  damping  to 
give  asymptotic  stability  (trajectories  spiral  to  an 
equilibrium  point  that  is  not  necessarily  the  origin) . 

Interestingly,  a  limit  cycle  exists  for  every  case  with 
negative  damping  (B<0) .   If  the  position  feedback  gain  K  is 
sufficiently  large  to  give  the  linear  system  a  complex  pair 
of  poles,  an  unstable,  unique  limit  cycle  exists,  whose 
amplitude  depends  on  K,  B  and  the  inertia  J  and  sliding 
friction  level  Lq.   This  "sliding"  type  of  friction  limit 
cycle  is  one  in  which  the  motor  or  other  object  being  moved 
never  "sticks";  the  reversing  torque  is  sufficient  when  it 
comes  to  rest  to  immediately  breakaway  in  the  opposite 
direction.   Therefore,  the  orbit  consists  of  two  parts,  each 
representing  motion,  with  two  switchings  per  cycle.   The 
initial  conditions  can  be  expressed  in  terms  of  the  system 
parameters  as 

(1.2)  xio  =  (Lc/K)  [(/x+l)/(M-l)],    X20  =  0 
where 

(1.3)  M  =  exp[-B7r/(2J;3)  ],  p    =    h    [4K/J  -  B^/J^]'"' 
Figure  1-4  shows  some  values  of  amplitude  (which  equals  x^^q) 
as  a  function  of  gain  K  for  the  case  J   =   Lq  =   1 ,    B  =  -1;  it 
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appears  that  a  limit  cycle  of  any  amplitude  can  be  found  by 
adjusting  gain  properly. 

The  unstable  limit  cycle  bounds  the  region  of 
asymptotic  stability:  trajectories  inside  spiral  in  to 
equilibrium  near  the  origin,  while  those  outside  go  to 
infinity  (Figure  1-5  shows  phase  plane  plot) .   This  example 
shows  why  it  is  dangerous  to  design  underdamped  systems,  and 
depend  on  the  sliding  friction  to  stabilize  the  system;  the 
friction  damping  is  amplitude-dependent,  and  at  some 
critical  amplitude,  is  insufficient  for  stability. 
Example  1-2;   Three-Dimensional  (3D)  System  with  Both 
Sliding  and  Sticking  Limit  Cycles 

The  second  example  system  to  be  considered  is  defined 
by  the  differential  equations 
(1.4)      dy/dt   =   Ay   +   b  Lp 
where 


(1.5) 


A   = 


0 
0 


1 

0 

-K, 


0 
1 
-B 


1       2 
with  b'  =  [  0  0-1  ],  Lp  =  friction  torque,  and  where  the 

state  variables  are  y  =  position,  y_  =  velocity,  and 

y   =  compensator  integrator.   Figure  1-6  shows  a  block 

diagram  of  this  system. 

Case  I:   One  Stable  Pole  and  One  Imaginary  Pole  Pair  (Only 

Sticking  Limit  Cycle  Exists) 

For  specific  results,  let  us  first  set  B  =  K^  ~  ^2  ~  "'" ' 

and  set  the  sliding  friction  torque  to  1.0  and  sticky 
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FIGURE  1-6.   BLOCK  DIAGRAM  OF  3D  SYSTEM  (EXAMPLE  2) 
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friction  (breakaway  torque)  to  1.2.   Simulation  of  this  case 
demonstrates  a  sticking  limit  cycle  with  an  initial 
condition  ofy'  =  [  1.0   0.2   0.0  ].   This  type  of  friction 
limit  cycle  repeats  a  four-part  cycle  of  sticking  (until 
torque  integrates  up  to  breakaway) ,  sliding  to  a  new 
position,  sticking  again,  and  sliding  back  to  the  start 
position.   A  visual  display  of  this  behavior  can  be  seen  in 
Figures  1-7  through  1-9,  which  show  projected  views  of  the 
closely  similar  limit  cycle  found  using  the  Case  II 
parameters. 

The  solution  of  the  piecewise-linear  equations  results 
in  a  valid  limit  cycle  solution  where  y^, '  =  [  1   0.2   0  ], 
T   (the  sliding  period)  =  it,    and  the  sticking  period 
T_  =  10,  yielding  an  overall  limit  cycle  period  of  approx. 
26.3  seconds,  matching  the  simulation  results. 
Case  II:   One  Stable  Pole  and  One  Unstable  Pole  Pair  (Both 
Sticking  and  Sliding  Limit  Cycles  Exist) 

If  the  same  example  system  is  used  with  B  =  0.9,  and 
K   =  K_  =  1  still,  the  closed-loop  eigenvalues  (of  the 
linear  portion  of  the  system)  move  into  the  right-half 
plane,  to  .026  ±  j  1.024,  while  the  other  pole  is  at  -.9524. 
As  might  be  expected,  there  is  still  a  sticking  limit  cycle 
solution  close  to  that  of  Case  I  (Figures  1-7  through  1-9) . 
Numerical  methods  give  a  solution  at  approximately 
Yq  =  [0.98   0.22   0  ],  and  approximately  the  same  period 
(confirmed  by  simulation) . 
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The  interesting  point  about  this  case,  however,  is  the 
similarity  to  Example  I-l,  the  2D  case;  that  is,  an  unstable 
complex  pole  pair  exists.   Since  the  nonlinear  damping 
effect  of  the  sliding  friction  is  largest  at  small 
amplitudes,  we  might  expect  trajectories  to  converge  for 
small  amplitudes  and  diverge  at  large  amplitudes  where  the 
unstable  linear  poles  overcome  the  nonlinear  damping  (as  in 
the  2D  case) .   Therefore,  this  case  would  have  both  a 
sticking  and  a  pure  sliding  friction  limit  cycle! 

Simulation  having  shown  a  limit  cycle  with  a 
half-period  T^  approximately  equal  to  3.1  seconds,  numerical 
solution  of  the  equations  for  the  symmetric,  sliding  case 
was  performed,  yielding  a  solution  at  T   =  3.1455  seconds, 
and  y   =  [-0.283   12.4   0.].   Figures  I-IO  and  I-ll  show 
some  views  of  the  two  3D  limit  cycles  for  this  case. 

It  is  expected  from  physical  intuition  and  simulation 
results  that  the  sticking  orbit  would  be  stable,  while  the 
sliding  orbit  would  be  unstable.   This  is  in  fact  the  case, 
with  the  sliding  orbit  exhibiting  saddle  point  behavior--one 
stable  and  one  unstable  mode  (real  eigenvalue  /  eigenvector 
pair) ,  while  the  sticking  orbit  has  stable  node  behavior-- 
two  real  stable  modes.   The  analysis  of  the  orbital 
stability  will  be  examined  in  Chapter  V,  where  Example  1-2, 
case  2  will  be  used  to  demonstrate  the  method  used  for  limit 
cycle  orbital  stability  calculations. 
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Organization  of  Dissertation 

The  main  body  of  this  dissertation  is  organized  as 
follows.   After  a  review  of  existing  results  on  this  problem 
and  related  areas,  the  standard  model  introduced  above  is 
developed  more  fully.   The  objective  in  the  model 
development  (in  Chapter  III)  is  to  provide  more  convenient 
forms  for  analysis  (a  control  canonical  form  and  a  normal 
form) ,  while  maintaining  the  generality  and  traceability  to 
actual  servomechanisms  of  the  original  model. 

Once  the  model  is  in  the  proper  form,  Chapter  IV 
develops  the  nonlinear  system  of  equations  that  represent 
the  necessary  and  sufficient  conditions  for  the  existence  of 
a  friction  limit  cycle.   Several  special  cases  are  examined 
that  are  somewhat  simplified  from  the  general  problem  and 
easier  to  solve.   For  example,  if  it  can  be  assumed  that  the 
limit  cycle  is  symmetric  about  the  origin  (that  is, 
half-wave  symmetric) ,  the  equation  system  can  be  simplified 
considerably. 

The  benefits  of  symmetry  in  the  solution  of  the  limit 
cycle  conditions  led  to  an  examination  of  the  conditions 
under  which  this  symmetry  existed.   Chapter  V  examines  this 
problem  and  describes  special  cases  for  which  symmetry 
exists.   Unfortunately,  it  is  not  known  if  friction  limit 
cycles  are  symmetric  in  general,  although  a  search  for 
counterexamples  was  unsuccessful. 

Using  the  results  of  the  previous  chapters,  a  complete 
stability  analysis  of  friction  limit  cycles  is  also  included 
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in  Chapter  V.   Using  the  eigenvalues  of  the  local  stability 
matrix  developed  in  this  analysis,  the  stability  of  limit 
cycles  previously  found  is  determined.   This  allows  the 
global  phase  portrait  of  the  system  to  be  pieced  together 
from  a  definition  of  all  the  equilibrium  points,  limit 
cycles,  and  their  stability  characteristics. 

The  results  presented  here  are  considered  only  a  start 
in  an  essentially  new  area:   use  of  piecewise  linear 
techniques  for  limit  cycle  analysis.   Many  possibilities  for 
further  research  exist,  both  in  the  completion  of  the 
solution  for  friction,  and  in  application  of  the  method  to 
other  piecewise  linear  nonlinearities.   As  pointed  out 
above,  the  majority  of  the  plant  nonlinearities  encountered 
in  servo  design  are  piecewise  linear,  including  gear 
backlash,  saturation,  quantization,  and  dead  zone.   The  last 
chapters  provide  a  summary  of  results  and  suggestions  for 
further  research  in  these  areas. 


CHAPTER  II 
REVIEW  OF  EXISTING  RESULTS 

An  Unpopular  Area  of  Study 
Surprisingly  little  work  has  been  done  in  this  area, 
for  various  reasons.   First,  the  nonlinearity  is  not  "nice," 
i.e.,  continuous  or  smooth,  which  prevents  the  application 
of  many  mathematical  tools.   Second,  when  modelled  with  two 
inputs  it  really  becomes  a  system  with  multiple 
nonlinearities,  for  which  standard  approaches  cannot  be 
used,  and  extremely  few  results  are  available. 

Finally,  and,  I  believe,  most  importantly,  it  is 
perceived  as  generally  benign  in  effect.   It  can  be  shown 
fairly  easily  (see  Thaler  and  Pastel,  1962,  for  example)  in 
the  2D  case  (i.e.,  two  state  variables),  that  if  the  linear 
portion  of  the  system  is  already  stable,  the  system  with 
pure  sliding  friction  added  is  also  stable.   It  seems 
probable  that  this  result  could  be  extended  to  arbitrary- 
dimensioned  systems,  based  on  the  following  energy  argument: 
the  linear  system  is  stable  (hence  contracting) ,  and  because 
the  (sliding)  friction  opposes  motion,  it  damps  the  system 
further,  so  the  map  is  more  contracting  than  without  it. 
Therefore,  friction  is  generally  held  to  be  benign  in  effect 
(it  increases  the  stability  of  a  system  by  increasing  the 
damping) . 
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For  these  and  possibly  other  reasons,  the  references 
available  in  this  area  were  few,  and  no  analysis  of  the 
general  multi-dimensional  friction  limit  cycle  problem  was 
found. 

In  spite  of  these  comments,  the  study  of  friction  limit 
cycles  has  practical  value  for  control  system  design. 
Although  the  effect  of  sliding  friction  is  to  increase 
damping,  it  has  become  fairly  standard  to  precompensate  for 
this  effect  in  a  design.   The  system  is  deliberately 
designed  to  be  underdamped,  so  that  it  is  not  too  well- 
damped  when  friction  is  added.   It  is  conjectured  that  this 
technique,  when  pushed  to  extremes,  results  in  large-signal 
instability,  since  the  damping  effect  varies  inversely  with 
amplitude. 

In  addition,  the  stiction  component  can  cause 
difficulties  also.   An  example  was  given  in  Chapter  I  of  a 
stable  (although  not  asymptotically  stable)  third-order 
linear  system  which  has  a  limit  cycle  when  nonlinear 
friction  (including  stiction  effects)  is  added.   A  slight 
variation  in  parameters  results  in  a  stable  linear  system 
which  destabilizes  with  the  addition  of  nonlinear  friction! 
Phase  Plane  (2D)  Analysis  of  Friction 
As  mentioned  above.  Thaler  and  Pastel  (1962),  in  their 
classic  text  on  nonlinear  systems,  completely  solved  the 
second-order  case  for  friction,  including  both  sliding  and 
static  components.   They  also  give  an  exact  criterion  for 
the  existence  of  limit  cycles  when  an  input  ramp  is  present. 
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They  showed  that  the  system  had  no  limit  cycles  (for  a 
system  with  positive  damping)  for  the  zero-input  case. 
Earlier  work  by  the  same  authors,  (Pastel  and  Thaler,  1960) 
actually  shows  a  stability  boundary  for  backlash  (as  a 
function  of  system  damping) ,  then  demonstrates  the 
stabilizing  effect  of  coulomb  friction. 

Analyses  of  the  2D  friction  problem  were  presented  in 
several  other  references  from  the  period  when  phase  plane 
analysis  was  an  active  area  of  research. 

Literature  on  Piecewise  Linear  Method 

The  approach  used  in  this  dissertation  to  give  exact 
solutions  for  limit  cycles  (if  they  exist)  has  been  known 
and  used  for  many  years.   In  the  relay  servomechanism 
literature,  it  is  known  as  the  piecewise  linear  method. 

The  method  was  used  primarily  for  relay  systems,  though 
it  is  generally  applicable  to  any  piecewise  linear 
nonlinearity .   Note  that  the  piecewise  linear  method  is 
essentially  the  same  as  phase  plane  analysis  for  2D  systems, 
so  there  is  overlap  between  the  phase  plane  references,  such 
as  those  cited  previously,  and  those  for  the  piecewise 
linear  method.   Unfortunately,  the  graphical  methods  of 
phase  plane  analysis  cannot  be  easily  extended  to  higher 
dimensions  (see  comment  below  on  the  work  of  Ku) . 

The  piecewise  linear  method  was  used  at  least  as  early 
as  1963  (by  Kovatch) ;  additional  studies  can  be  found  in 
references  such  as  Kovatch  (1964)  for  two  nonlinearities, 
O'Donnell  (1964)  for  time-optimal  switching,  Marstrander 
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(1969)  for  backlash,  Negoescu  and  Sebastion  (1971) ,  Langill 

(1965),  and  Urabe  (1967).   Weischedel  (1973)  applied  the 

method  to  on-off  control  systems.   Ku  (1958)  made 

significant  and  early  contributions,  including  attempts  to 

extend  the  phase  plane  method  to  higher  dimensions  by 

various  graphical  projections. 

The  Russian  controls  literature  calls  this  the  Andronov 

point  transformation  method;  a  reference  was  found  in  a 

translation  of  material  developed  in  1956  by  E.P.  Popov 

(translated  1962) .   Although  significant  development  of  this 

area  has  been  performed,  especially  in  the  Russian  journal 

Automatika.  an  application  to  the  friction  problem  has  not 

been  found. 

Extensions  to  Multi-Dimensional  Case 
(Approximate  Analvses  by  Describing  Function) 

Analysis  using  approximate  techniques  such  as  the 
describing  function  go  back  at  least  as  far  as  the  piecewise 
linear  method,  or  even  farther  in  the  Russian  journals 
(Popov,  1956).   Four  papers  are  listed  by  Gibson  (1963), 
that  he  refers  to  as  independent  developments  of  the 
describing  function  method;  all  are  from  the  late  1940s! 

However,  the  two  input  difficulty  that  limits  the 
applicability  of  describing  functions  to  friction  is  usually 
dealt  with  by  ignoring  the  static  friction.   No  references 
have  been  found  that  discuss  the  general  case  for  an  exact 
friction  model.   No  other  references  were  found  that 
attempted  the  multi-dimensional  friction  problem  with  other 
methods. 
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General  Results  on  Periodic  Orbits 

In  the  mathematics  literature  on  the  theory  of  periodic 
solutions  of  ordinary  differential  equations,  the  Poincare 
map  (Hirsch  and  Smale,  1974)  is  frequently  used;  this  map 
represents  the  change  in  system  state  after  one  cycle  or 
period.   By  searching  for  fixed  points  of  the  map  (which 
represent  a  return  to  the  original  state  after  one  cycle) ,  a 
periodic  solution  can  be  found.   This  theory  was  used  in 
generating  the  stability  results  to  be  stated  in  Chapter  V, 
where  the  orbital  stability  of  friction  limit  cycles  is 
examined.   The  Chapter  V  results  are,  therefore,  original 
only  in  their  application  to  friction  limit  cycles,  not  in 
the  theory  of  orbital  stability. 

There  are  many  results  in  this  literature  (although  the 
Poincare  theory  itself  is  restricted  to  planar  systems) . 
Pliss  (1966)  stated  some  results  on  the  existence  and 
uniqueness  of  these  fixed  points,  although  they  could  not 
apply  to  friction  as  they  were  restricted  to  continuous 
functions,  or  non-autonomous  systems.   Hale  (1969)  discussed 
fixed-point  theorems  for  2D  systems  and  n-dimensional 
systems  on  a  torus;  in  addition  he  stated  the  orbital 
stability  result  (that  all  eigenvalues  of  the  map  except  one 
must  be  less  than  one  in  magnitude) .   Hayashi  (1964)  refers 
to  this  approach  as  the  topological  approach  to  periodic 
solutions,  since  the  theorist  is  concerned  with  the 
existence  and  topological  properties  of  integral  curves  in 
state-space.   Finally,  Hirsch  and  Smale  (1974)  provide  an 
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excellent  discussion  on  an  introductory  level  of  the  theory 
of  flow  maps;  this  reference  was  invaluable  in  the  stability 
analysis  of  Chapter  V. 

My  survey  of  this  literature  was  thorough,  but  not 
exhaustive;  no  applications  to  this  nonlinearity  were 
encountered. 

Summary  of  Literature  Review 

In  the  20  years  since  the  heyday  of  phase  plane  and 
piecewise  linear  methods,  it  appears  that  almost  no  new 
results  have  been  discovered  in  this  area.   However,  the 
study  of  friction  limit  cycles  has  significant  practical 
value  for  designers  of  servomechanisms,  since  a  servo 
designed  using  standard  techniques  may  exhibit  instabilities 
in  the  presence  of  friction. 

The  literature  survey  was  extensive  in  the  area  of 
controls,  and  moderate  in  the  mathematics  field.   In  the 
former,  IEEE  indices  were  searched,  including  25  years  of 
the  Transactions  on  Automatic  Control.   Twenty-five  years  of 
the  Russian  journal  Automation  and  Remote  Control  were 
examined.   The  journal  Automatica  was  also  reviewed. 
Although  no  journals  in  the  field  of  pure  mathematics  were 
checked,  many  texts  in  the  areas  of  stability,  limit  cycles, 
and  dynamical  systems  were  examined. 

The  main  contribution  of  this  research  is  the 
application  of  known  (although  relatively  obscure  in 
controls  literature)  techniques  such  as  the  piecewise  linear 
method  and  Poincare  flow  maps  to  the  problem  of  nonlinear 


29 
friction.   No  reference  was  found  that  considers  the  multi- 
dimensional friction  limit  cycle  problem  using  an  exact 
method  such  as  the  piecewise  linear  technique.   All  the 
references  found  either  analyze  the  2D  case  for  friction,  or 
derive  approximate  results  by  using  a  simpler  friction 
model.   Practical  experience  with  servos  indicates  that  it 
is  unacceptable  to  ignore  either  friction  component. 

The  development  presented  here,  therefore,  appears  to 
be  the  first  to  attack  the  multi-dimensional  friction  limit 
cycle  problem  with  sufficient  accuracy  for  practical 
applications.   In  addition,  the  results  on  limit  cycle 
symmetry  for  piecewise  linear  systems,  and  the  derivation  of 
limit  cycle  orbital  stability  by  the  same  method,  are 
original,  to  the  author's  knowledge. 


CHAPTER  III 
DEVELOPMENT  OF  STANDARD  MODEL 

This  chapter  is  concerned  with  setting  up  the  model  of 
the  system  to  be  studied,  and  converting  it  to  a  convenient 
form  for  analysis.   The  original  system  model  is  in  the  form 
of  physical  state  variables,  where  the  states  represent 
actual  physical  quantities  in  the  servomechanism.   The  model 
form  to  be  used  for  most  results  is  a  control  canonical 
representation,  which  is  developed  from  the  physical  states 
by  a  linear  transformation.   This  form  has  advantages  in 
clarifying  subsequent  developments  of  the  theory,  including 
the  role  of  system  zeroes  and  characterizing  limit  cycle 
symmetry.   This  chapter  also  demonstrates  how  the  piecewise 
linear  friction  torque  input  can  be  removed  by  a 
translation,  which  leads  to  a  simplification  of  the 
existence  formulas  of  Chapter  IV. 

In  addition  to  the  derivation  of  the  state  forms,  this 
chapter  provides  a  description  of  the  various  assumptions 
made  in  the  analysis.   In  other  words,  the  generality  of  the 
analysis  and  the  solutions  are  defined  by  examining  the 
assumptions  made.   For  example,  although  the  original 
motor/load  model  included  no  compliance  or  resonance 
effects,  this  restriction  can  be  removed,  allowing  the 
results  to  be  more  generally  applicable. 
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System  Model  with  Physical  State  Variables 

The  two  primary  physical  variables  are  angular  velocity 
(required  as  input  to  the  sliding  mode  friction)  and  its 
integral,  angular  position.   The  source  of  this  model  is  a 
rotating  servomechanism.   The  theory  is  not  restricted  to 
rotational  friction  problems,  however,  since  angular 
variables  (torque,  angular  velocity,  etc.)  could  be  replaced 
by  linear  ones  (force,  velocity,  etc.)  without  affecting  the 
nature  of  the  problem  (this  point  is  discussed  further  at 
the  end  of  this  chapter) .   The  model  was  illustrated  in 
block  diagram  form  in  Figure  I-l;  the  controller  transfer 
function  G(s)  is  replaced  by  state  equations  in  the  model 
below. 

The  remaining  n-2  state  variables,  where  n  is  the  order 

of  the  system,  are  controller  states.   That  is,  they  form 

the  model  of  the  servo  loop  compensation,  amplifiers,  motor 

armature  effects,  and  other  dynamic  portions  of  the  system. 

Thus  we  may  set  up  a  model  of  the  following  form: 

(3.1)      Y=AY+a-y,+a_y 
^    '      -^c    c-^c   — cl-'n-l   -c2-'n 

^n-1  "  ^n 

J  y  =  -K^y   ^  -  By   +  b'y  +  f (v) 
■'n     l-'n-l    -'n   -  -^c     ^^' 

where  the  y.  variables  are  the  physical  states  and  make  up 
the  vector  y,  the  dot  (')  indicates  the  time  derivative 
d/dt,  prime  (')  indicates  transpose,  f(.)  is  the  friction 
torque,  J  is  the  lumped  system  inertia,  and  B  is  the 
damping.   The  n-2  dimensional  linear  subsystem  in  the  n-2 
states  in  vector  y   is  formed  by  the  system  matrix  A 
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(n-2  X  n-2) ,  the  input  vectors  a    and  a    (n-2  x  1) ,  and 

the  output  vector  b  (n-2  x  1)  that  feeds  the  output  of  the 

controller  into  the  torque  equation. 

The  only  assumptions  made  in  the  construction  of  this 

model  are  that  y   .  and  y   are  angular  position  and 

n-1       n        ^      ^ 

velocity,  respectively.   These  assumptions  completely  define 
the  first  equation,  and  require  that  the  nonlinear  friction 
torque  term  f/J  appears  in  the  y  differential  equation 
(representing  angular  acceleration  due  to  the  sum  of  torques 
applied) .   The  rest  of  the  system  is  linear;  the  friction 
term  does  not  appear  in  any  other  derivative. 

Note  that  there  is  an  assum.ption  made  here  that 
restricts  the  generality  of  the  model:  that  the  friction 
appears  in  only  one  equation.   Although  examples  can  be 
constructed  that  do  not  obey  this  restriction  (such  as  those 
involving  differentiators,  full  state  feedback,  torque 
measurement,  or  a  foundation  model) ,  this  model  is  felt  to 
be  quite  general,  and  covers  essentially  all  servo  problems 
where  friction  limit  cycle  information  is  of  interest.   This 
point  will  be  discussed  in  more  detail  later  in  this 
chapter. 

Modelling  Friction  By  a  Nonlinear,  Two-Input  Function 
The  nonlinearity  investigated  in  this  analysis  is  a 
model  for  mechanical  friction  that  depends  on  both  the 
relative  velocity  of  the  surfaces  and  the  driving  torque  of 
the  moving  element.   The  level  of  friction  opposing  the 
motion  depends  on  velocity  (or  more  precisely,  the  direction 
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of  the  velocity)  when  the  element  is  moving.   In  the  case 
where  the  element  is  stationary,  on  the  other  hand,  the 
static  friction  torque  holding  the  system  motionless  must  be 
set  equal  to  the  input  driving  torque  that  is  attempting  to 
move  the  element.   Each  of  these  cases  is  based  on  the 
models  of  Coulomb  friction  described  in  any  basic  statics 
textbook.   Figure  1-2  in  Chapter  I  illustrates  the  friction 
model . 

Although  an  analysis  can  be  performed  with  only  one 
element  modelled,  thereby  using  a  single-input  nonlinearity 
and  simplifying  the  analysis,  this  would  lead  to  results 
that  would  not  apply  in  any  physical  system,  and  could  be 
misleading  or  erroneous.   Cases  can  be  constructed  where  the 
modelling  of  sliding  friction  alone  predicts  a  lim.it  cycle, 
due  to  its  highly  nonlinear,  discontinuous  characteristic. 
In  a  practical  system,  however,  no  limit  cycle  would  exist. 

For  example,  the  two-dimensional  friction  limit  cycle 
is  completely  solved  in  Thaler  and  Pastel  (1962) ,  as 
described  in  Chapter  I.   Neglecting  the  static  friction 
component  in  this  analysis  leads  to  the  conclusion  that  a 
limit  cycle  of  very  small  amplitude  can  be  obtained. 
However,  in  reality,  such  a  small  initial  condition  will 
stick  at  the  initial  point  due  to  the  static  friction 
torque. 

On  the  other  hand,  the  sliding  component  of  friction 
cannot  be  ignored  either.   As  a  qualitative  example  of  a 
case  demonstrating  this,  consider  a  system  with  no  damping. 
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which  would  limit  cycle  if  static  friction  alone  were 
present.   The  system  starts  at  an  initial  position  with 
nonzero  error,  integrates  up  in  the  controller  until  it 
breaks  free  from  the  static  friction,  and  overshoots  the 
desired  position  to  start  the  cycle  all  over  again.   The 
inclusion  of  sliding  friction  in  the  model,  however,  might 
provide  enough  damping  to  the  system  to  cause  it  to  spiral 
in  to  the  equilibrium  point  at  the  origin.   Another 
possibility  indicated  by  the  behavior  of  typical  friction 
limit  cycles  found  is  that  a  system  that  is  actually  stable 
for  small  signals  (bounded  region  of  stability)  might  appear 
unstable  if  analyzed  without  sliding  friction. 

These  cases  indicate  that  accurately  predicting  limit 
cycles  for  physical  servomechanisms  requires  a  complete 
model  of  friction,  as  is  used  in  the  analysis  contained  in 
this  paper. 

The  Physical  State  Model  in  Piecewise  Linear  Regions 

As  described  in  the  introduction,  application  of  the 
piecewise  linear  method  to  the  system  model  results  in  five 
regions  of  state  space,  with  a  linear  system  for  each 
region.   The  form  of  the  (now  linear)  model  in  each  region 
is  described,  so  that  the  trajectories  in  each  region  can  be 
described. 

Region  I:  Velocity  <  0  (y  <0) .   This  n-dimensional  region  is 
described  by  the  same  set  of  differential  equations  as  in 
(3.1),  except  that  the  nonlinear  friction  term  f(Y)  is 
replaced  by  a  constant  input  L  .   This  follows  from  the  fact 


35 
that  the  friction  in  sliding  mode  is  constant  (independent 
of  velocity)  as  long  as  the  system  is  in  motion  in  one 
direction.   From  basic  mechanics,  this  constant  friction  is 
a  function  of  the  normal  force,  with  the  proportionality 
constant  being  the  coefficient  of  sliding  friction.   The 
torque  L   is  positive,  since  physical  friction  always  acts 
to  oppose  the  direction  of  motion;  velocity  is  negative,  so 
the  acceleration  torque  due  to  friction  is  positive. 
Therefore,  the  model  in  region  I  is  a  linear  system  driven 
by  a  constant  input. 

Region  II;  Velocity  >  0  (y  >0) .   This  region  model  is 
exactly  the  same  as  in  region  I,  except  that  f(y)  is 
replaced  by  -L  . 

Region  III;  Velocity  =  0  (y  =0)  and  Acceleration  Torque  < 
Breakaway  Torque.   In  this  case,  the  system  is  in  a  static 
condition,  and  the  model  for  the  static  component  of 
friction  applies  (also  known  as  stiction) .   The  friction 
opposes  the  applied  acceleration  torque  with  an  equal  torque 
so  that  the  net  acceleration  torque  is  zero,  up  to  a  maximum 
amount  of  static  friction.   The  maximum  amount  of  friction 
(the  breakaway  torque)  is  determined  by  the  coefficient  of 
static  friction  and  the  normal  force.   Once  the  acceleration 
torque  becomes  greater  than  the  maximum  stiction,  the  net 
torque  becomes  nonzero,  and  the  system  begins  to  move 
(enters  region  I  or  II) . 

Note  that  the  acceleration  torque  used  to  define  this 
region,  and  which  is  matched  by  an  opposing  friction  force. 
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consists  of  every  term  in  the  torque  equation  (the 

derivative  of  y  )  except  the  friction  f.   Therefore,  the 

linear  terms  in  the  torque  equation  define  the  boundaries  of 

this  region.   The  region  is  contained  in  the 

(n-1) -dimensional  hyperplane  y_,=0- 

Since  the  net  torque  is  zero  and  the  velocity  is  zero, 

the  first  two  differential  equations  in  (3.1)  drop  out  (the 

state  y   ,  remains  constant  and  y   remains  zero  over 
■'n-l  -'n 

trajectories  in  this  region) ,  leaving  the  n-2  dimensional 

controller  system  as  the  only  equations  required  to  define 

the  trajectory  in  region  III.   Of  course,  input  vector  a 

drops  out  (since  yn=0) /  while  a    and  b  are  used  to  model 

the  effect  of  the  constant  angular  position  on  the 

controller  and  the  controller  output  torque  (to  determine 

the  breakaway  condition) . 

Region  IV:  Velocity  =  0  (y  =0)  and  Acceleration  Torque  > 

n 

(+) Breakaway  Torque.   This  region  is  also  contained  in  the 
(n-1) -dimensional  hyperplane  Yx^=0 ,    but  the  driving 
acceleration  is  sufficient  at  any  point  to  overcome  the 
maximum  static  friction  opposing  the  impending  motion.   Thus 
a  system  that  comes  to  rest  with  sufficient  torque 
(trajectory  enters  region  IV  rather  than  III)  immediately 
breaks  into  motion  again,  entering  region  II  (since 
acceleration  torque  is  positive,  velocity  increases  from 
zero  and  trajectory  must  enter  region  II). 

For  this  reason,  a  trajectory  can  be  said  to  cross  this 
region  (in  the  sense  of  crossing  a  plane,  for  example),  but 
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is  only  in  the  region  for  an  instant.   Since  the  integral  of 
a  finite  quantity  over  a  set  of  measure  zero  is  zero,  the 
applied  torque  and  the  differential  equations  are  irrelevant 
to  the  system  behavior.   The  trajectory  leaves  the  region  at 
the  exact  point  it  entered.   It  is  therefore  unnecessary  to 
examine  the  form  of  the  system  equations  in  this  region. 
Region  V:  Velocity  =  0  (v, '^Q)  and  Acceleration  Torque  < 
(-)  Breakaway  Torque.   The  region  V  model  is  identical  to 
the  region  IV  model,  except  for  the  breakaway  conditions. 

The  five  regions  cover  the  state  space,  including  the 
equilibrium  point  (actually,  line  segment)  at  the  origin. 
The  trajectories  in  each  of  the  individual  regions  can  be 
described  by  the  solution  to  the  linear  system  of 
differential  equations: 
(3.2)      Y(t)  =  F(t,tQ)  y(tQ)  +  Git,t^)    L 

F(t,tQ)  =  exp[A(t-tQ)] 
where  L  is  the  (constant)  input  in  that  region,  A  is  the 
system  matrix  for  that  region,  and  the  F  and  G  matrices  vary 
from  region  to  region.   The  initial  conditions  Y(t  )  are 
determined  by  the  commutation  conditions  at  the  boundary  of 
the  region  at  which  the  trajectory  enters. 

An  algebraic  trick,  to  be  presented  below,  can  be  used 
with  this  form  to  remove  the  input  term  (GL)  in  equation 
(3.2).   The  system  then  behaves  as  an  autonomous  system, 
with  considerable  simplification  of  the  defining  equations. 
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Conversion  to  Control  Canonical  Form 
The  most  convenient  form  for  analysis  is  the  control 
canonical  form,  because  the  role  of  system  poles  and  zeroes 
in  the  existence  conditions  in  Chapter  IV  will  be  clarified 
through  the  use  of  this  representation.   In  addition,  a 
result  in  Chapter  V  on  symmetry  is  proved  using  this  form. 
Although  this  section  constitutes  a  digression  from  the  main 
line  of  development,  it  is  justified  by  the  usefulness  of 
the  control  form. 

The  conversion  to  control  form  is  accomplished  by 
representing  the  controller  (transfer  function  G(s)  in 
Figure  I-l)  in  control  form  and  combining  this 
(n-2) -dimensional  system  with  the  equations  for  position  and 
velocity.   When  completed,  the  system  model  will  have  the 
following  form: 
(3.3)      x^  =  x^ 

•      •      •      • 

X   -  =  X 

n-1    n 

x_  =  "Q:  X   -  a^x^  -  ...   -  a  X   +  f(x) 
n      0  112  n  n     ^—' 

The  coefficients  in  the  differential  equation  for  x   are  the 
coefficients  of  the  closed  loop  polynomial  of  the  linear 
system. 

The  procedure  used  to  convert  the  system  to  control 

canonical  form  is  as  follows.   The  physical  state  model  is 

set  up  so  that  the  controller  (the  (n-2) -dimensional 

subsystem  represented  by  state  variables  y  )  is  in  control 

c 
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canonical  form  already.   This  can  be  done  for  any  controller 
representation  that  is  controllable  (if  the  original 
representation  for  the  controller  was  a  transfer  function,  a 
controllable  representation  can  always  be  found) . 

Next,  the  model  is  transformed  so  the  velocity  feedback 
into  the  controller  subsystem  is  zero  (i.e.,   vector  a    in 
(3.1)  is  zero).   This  can  be  done  without  loss  of  generality 
in  the  model,  as  shown  by  the  following  argument:   Suppose 
a    is  not  zero  in  the  original  system.    We  can  apply  a 
similarity  transformation 
(3.4)      Y  =  T^  y(l) 
where  yC^)  is  the  new  state  variable  vector,  and 


T   = 
a 


10  0 
0  10 


0  a 


c2 


0  a 


c2 


0 
0 


00.   .   .   0   la^"^0 

C2 

0  0  ...   0   0    1     0 
0   .    .    .     0    0     1 
Note  this  is  equivalent  to  replacing  all  (n-2)  controller 

states  Yc   ^Y    (Yc^'''^  +  ^c2  Vn-l)  •   •^  straight  forward 
calculation  (plug  (3.4)  into  (3.1))  shows  that  the  resulting 
system  has  no  feedback  of  velocity  into  the  controller. 

Since  a  similarity  transformation  (or  the  equivalent 
change  of  variables)  does  not  affect  the  basic  system 
behavior,  but  merely  its  representational  model,  this 
alteration  of  the  model  is  valid.   Physically,  this 
operation  has  removed  a  redundant  parameter  in  the  original 
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model,  since  the  velocity  feedback  in  the  original  model  can 
be  achieved  by  adjusting  the  damping  parameter,  B,  or 
adjusting  the  position  feedback  gain,  K^ .   Since  the  models 
are  equivalent,  the  superscript  on  the  y  state  variables  can 
then  be  dropped  in  subsequent  equations. 

With  the  controller  in  canonical  form,  and  the  velocity 
feedback  to  the  controller  eliminated,  the  following 
intermediate  form  is  obtained: 


(3.5) 


y   =  A  y   +  a  ,y 
•^c     c-^c    -cl-^n-1 


^n-1  =  yn 


^  ^n  =  -^l^n-l  -  B^n  +  ^'^^c  +  ^^^^ 
In  order  to  complete  the  conversion  of  the  model  to 

control  form,  two  similarity  transformations  are  then 

applied  to  this  intermediate  model  in  succession.   As 

previously  explained,  similarity  tranf ormations  do  not 

change  the  basic  behavior,  hence  are  valid  alterations  of 

the  model  equation.   The  transformations  are: 

(3.6a) 

(3.6b) 

where  the  nxn  matrix  A  is  obtained  from  the  model  in  (3.5), 
and 


A^  =  T^  ■*"  A  T^ 
^2  =  ^2~^   ^1  ^2 


(3.7a) 


^1  = 


1       0   . 
0   10 


0   0 
0 


0        .        .        . 

0       10 

0 

^  ^1  •   • 

•      "n-S    1 

0 

0       .           .           . 

0       0 

1 
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(3.7b)      T^  = 


1   0   . 
0   10 


0   0 
0 


0  ^0  ^,  . 


^n-3  1 


and  where  the  p . ' s   are  the  coefficients  of  the 
characteristic  polynomial  of  the  controller  alone  (i.e.,  the 
characteristic  polynomial  of  matrix  A   in  equation  3.5). 
The  transformed  system  is  then  in  control  canonical  form. 
Note  that  the  friction  input  is  unchanged  by  these 
transformations,  and  is  still  applied  in  the  y   differential 
equation.   The  composite  transformation  is  required  for 
later  derivations,  as  it  defines  the  relationship  between 
the  physical  states  and  the  control  form  states: 
(3.8)      y  =  T^T^x 
where 

10...      GO 

0   10...       0 


T    =  T  T   = 
12    ■'l'-2 


0   0 


^0  ^1  • 
°   ^0   ^1 


^.- 


n-3 


^n-3  1 


Simplification  of  DE  Solution  Via  Coordinate  Translation 

Once  the  system  is  in  control  form,  the  following 

translation  of  the  state  coordinates  can  be  used  to  convert 

it  to  an  autonomous  system  (within  each  region) .   Since  the 


42 


friction  is  a  constant  in  region  I  or  II  (L  or  -L  ) ,  and 
since  none  of  the  differential  equations  except  that  for  x 
involves  x  ,  the  x   coordinate  can  be  translated  by 

(3.9)  x^^  =  x^  -  f/a^ 

where  the  subscript  "a"  indicates  autonomous  coordinates, 

and  f  is  the  (constant)  friction  level  in  that  region  (I  or 

II) .   Substitution  of  this  translation  shows  that  the  x 

n 

equation  still  gives  the  same  result,  the  system  matrix  is 
unchanged,  and  the  input  has  disappeared.   The  autonomous 
system  model  is  then 

(3.10)  x^^  =  x^ 

*2  =  ^3 


x   ,  =  x 
n-1    n 

x^  =  -q:_Xt  -   a^x^    -    .  .  .      -  a  X 
n     0  la    1  2  n  n 

Therefore,  by  performing  a  translation  on  the  initial 
condition  and  on  the  final  state  of  the  trajectory  in  the 
region,  the  trajectory  can  now  be  described  by  the  simpler 
form 

(3.11)     x^(t)  =  F(t,tQ)  K^(t^) 
F(t,tQ)  =  exp[A(t-tQ)] 
instead  of  equation  (3.2),  where  the  matrix  A  is  in  control 
form  (i.e.,  the  matrix  A  produced  by  transformation  3.6). 

Of  course,  the  control  form  representation  is  not 
required  to  perform  this  translation.   The  equivalent 
translation  in  terms  of  the  physical  state  variables  can  be 
defined  using  the  transformation  matrix  as 
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(3.12)  ly  =  T^2  i 

Ya  =  i^  -  iy 
where  1'  =  [  ^/ccq      0  ...  0  ]  is  the  translation  of  the 
control  form  states  (equation  (3.9))  and  y^  is  the 
autonomous  physical  state  variables.   The  resulting 
differential  equation  for  y  would  be  of  the  form 

(3.13)  Ya  =  A  ya 

(with  A  unchanged) .   A  slightly  different  form  of  the 
existence  conditions  of  Chapter  IV  would  then  result.   The 
control  form,  has  certain  advantages,  however,  as  will  be 
seen. 
Summary  of  Assumptions  and  Discussion  of  Generality  of  Model 

The  following  assumptions  and  conditions  on  the  system 
are  made  during  the  analysis  of  friction  limit  cycles  in 
this  dissertation: 

(Al)   The  system  can  be  modelled  by  a  model  of  the  form 
in  Figure  I-l,  and  equations  (3.1). 

(A2)   The  velocity  feedback  into  the  controller  is  zero 
(without  loss  of  generality,  as  shown  by  argument  given 
previously  in  this  chapter) . 

(A3)   The  controller  portion  has  a  controllable 
representation  (this  is  assured  if  the  controller  is 
representable  by  a  strictly  proper  transfer  function;  if 
transfer  function  is  proper  the  direct  feed-through  term  can 
be  divided  out  and  lumped  with  the  direct  position  feedback, 
so  this  case  is  also  allowed) . 
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(A4)   Friction  appears  only  in  the  differential 
equation  for  velocity  (the  acceleration  torque  equation) . 
(A5)   The  closed-loop,  linear  portion  of  the  system 
should  not  have  a  pole  at  the  origin. 

(A6)   The  limit  cycles  to  be  analyzed  are  all  simple, 
that  is,  they  traverse  regions  I  and  II  once  before 
returning  to  the  initial  point,  and  therefore  have  only  four 
switchings  per  cycle  at  most. 

Note  that  many  of  the  assumptions  were  made  for 
convenience  in  the  subsequent  development,  so  that  a 
specific  model  could  be  used  for  concreteness ,   Most  can  be 
removed  by  a  transformation  of  the  original  problem  as 
discussed  below,  or  by  minor  modifications  of  the  theory, 
and  so  involve  no  loss  of  generality. 

In  fact,  the  only  restriction  to  the  system  model  that 
cannot  be  removed  is  that  the  limit  cycle  is  simple 
(assumption  A6) .   The  analysis  could  be  extended,  as 
discussed  below,  to  remove  even  this  restriction.   Thus  the 
theory  is  quite  general. 

Several  of  these  assumptions  require  further 
discussion.   Assumption  (A4)  is  required  to  perform  the 
analysis  using  the  piecewise  linear  method  in  Chapter  IV, 
since  the  analysis  uses  the  control  canonical  form  heavily. 
If  the  friction  input  appears  in  more  than  one  differential 
equation,  the  trick  used  there  to  eliminate  the  friction  by 
a  coordinate  translation  is  invalid. 
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Examples  could  be  generated  that  violate  this 
assumption.   If  acceleration  torque  was  directly  measured, 
the  friction  torque  could  be  inferred  and  used  as  an  input 
to  the  controller.   The  same  comment  holds  for  systems  with 
differentiators,  since  the  derivative  of  the  velocity 
involves  the  friction.   Friction  also  causes  a  reaction 
torque  on  the  mounting  of  the  motor  or  other  physical 
device,  so  if  a  foundation  or  mount  model  is  included,  the 
assumption  would  be  violated. 

This  condition  does  not  restrict  the  generality  of  the 
method,  however.   Systems  which  violate  this  assumption  can 
still  be  transformed  to  control  form,  the  method  applied, 
and  existence  conditions  checked.   A  minor  modification  of 
the  Chapter  IV  conditions  (involving  the  §_   vector  only) 
would  be  required. 

Assumption  (A5)  is  not  necessary,  except  for  steps 
involving  taking  the  inverse  of  certain  matrix  quantities  in 
Chapter  IV.   However,  the  assumption  results  in  no  loss  of 
generality,  for  the  following  reason:   a  system  with  a  zero 
eigenvalue  has  a  free  integrator  so  that,  with  the  proper 
selection  of  state  variables,  one  state  does  not  feed  back 
into  the  system.   This  can  be  seen  by  examining  the  control 
form  of  the  system  model;  when  an  eigenvalue  is  zero,  the 
coefficient  ag  of  the  characteristic  polynomial  is  zero,  and 
the  first  column  of  the  system  matrix  is  zero  (see  equation 
3.3).   Therefore  the  magnitude  of  the  state  x^  has  no  effect 
on  the  behavior  of  the  system.   Any  limit  cycle  existing  in 
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the  system  will  thus  also  exist  in  the  (n-1) -dimensional 
system  formed  by  deleting  the  state  x^  (although  not 
necessarily  vice  versa) .   Therefore,  any  system  that  fails 
this  assumption  can  be  analyzed  by  deleting  free  integrator 
states,  then  looking  for  limit  cycles;  any  limit  cycles 
existing  in  the  original  system  would  be  found  by  this 
method.   Note,  however,  that  solutions  may  be  introduced 
that  do  not  exist  in  the  original.   For  example,  an 
asymmetric  limit  cycle  in  the  reduced  system  might  not  exist 
in  the  original  since  the  free  integrator  output  might  not 
be  periodic  (see  Chapter  V) . 

Assumption  (A6)  is  also  necessary  to  allow  the 
derivation  of  Chapter  IV  to  proceed.   The  piecewise  linear 
method  requires  the  overall  shape  of  the  limit  cycle  to  be 
known,  including  which  regions  are  traversed  and  the  order 
of  traversal.   This  could  be  considered  a  limitation  of  the 
method,  since  the  analysis  cannot  be  performed  without 
knowing  the  shape  of  limit  cycles  that  are  to  be  found  by 
the  analysis! 

Note,  however,  that  an  even  more  restrictive  assumption 
is  made  in  the  classical  describing  function  analysis,  a 
tool  that  has  proven  quite  useful  for  decades.   That 
analysis  assumes  a  sinusoidal  input  to  the  nonlinearity , 
therefore  implicitly  assuming  a  half-wave  symmetric  type  of 
solution  with  two  switchings  per  cycle.   It  is  therefore 
felt  that  the  assumption  is  reasonable.   Future  research 
could  remove  this  assumption  by  listing  all  possible 


47 
solution  shapes  (and  order  and  number  of  traversals  of 
regions  during  each  cycle) ,  and  applying  the  method  to  each 
possibility.   The  difficulty  lies  in  the  large  number  of 
potential  solutions  to  be  checked. 

If  it  could  be  proved  that  all  limit  cycles  are  simple, 
the  difficulty  would  disappear.   Without  such  a  proof,  the 
results  of  this  thesis  apply  only  to  simple  limit  cycles. 
For  example,  the  necessary  conditions  for  the  existence  of  a 
limit  cycle  found  in  the  next  chapter  are  necessary  only  for 
simple  limit  cycles.   Therefore,  the  lack  of  a  solution  for 
these  conditions  for  a  particular  system  only  guarantees 
that  a  simple  limit  cycle  cannot  exist  for  that  system. 
More  complicated  behaviors  are  not  ruled  out. 

Finally,  some  comments  on  the  generality  of  the  system 
model  assumed  are  in  order  (assumption  (Al) ) .   The  model  was 
inspired  by  the  case  of  a  rotational  servomechanism,  such  as 
a  DC  motor.   Although  this  is  an  important  application,  the 
theory  developed  here  can  be  applied  to  many  other  systems, 
as  long  as  they  can  be  put  in  the  form  specified.   For 
example,  a  linear  drive  system  has  the  same  characteristics 
of  inertia,  damping,  friction,  and  acceleration  force 
(instead  of  torque) ,  so  the  theory  can  be  directly  applied. 
Systems  with  gear-trains  and  loads  can  be  modelled  by 
lumping  the  load  and  gear  inertia  into  the  motor,  as  long  as 
the  system  is  rigid  and  closely  coupled  (gear  backlash  can 
be  ignored) .   Although  the  model  is  for  a  continuous-time 
system,  the  use  of  discrete-time  (i.e.,  digital)  controllers 
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can  be  accomodated  by  using  their  continous-time 
equivalents. 

Many  systems  of  practical  interest  have  mechanical 
flexure,  modelled  by  resonances.   If  they  can  be  lumped  into 
the  model  of  the  controller,  the  theory  can  still  be  applied 
as  is.   However,  the  case  of  two  masses,  separated  by  a 
spring,  where  both  masses  experience  nonlinear  friction, 
does  not  fit  into  the  present  model.   Additional  piecewise 
linear  regions  must  be  specified,  with  a  great  increase  in 
the  complexity  of  the  analysis,  unless  the  limit  cycle  is 
known  to  be  low  in  frequency  so  the  system  could  be 
considered  to  behave  as  a  rigid  body. 

In  summary,  certain  assumptions  are  made  to  allow  a 
specific,  concrete  model  to  be  used  in  the  development. 
Examination  of  the  assumptions  shows  that  most  are  for 
convenience  only,  and  the  theory  applies  to  a  quite  general 
class  of  problems. 

Given  these  restrictions  and  assumptions  on  the  problem 
to  be  studied,  and  the  development  of  a  convenient 
representation,  the  piecewise  linear  method  can  be  applied 
to  the  problem  of  friction  limit  cycles.   The  next  chapter 
derives  exact  conditions  for  the  existence  of  these  limit 
cycles,  and,  once  a  solution  is  found,  provides  exact 
information  on  the  limit  cycle  trajectory. 


CHAPTER  IV 
EXISTENCE  CONDITIONS  FOR  LIMIT  CYCLES 

This  chapter  applies  the  piecewise  linear  method  to  the 
problem  of  friction  limit  cycles,  using  the  model  developed 
in  the  previous  chapter.   The  method  provides  a  set  of 
nonlinear  algebraic  equations  for  which  a  solution  must 
exist  in  order  to  have  a  limit  cycle.  These  represent 
necessary  conditions  for  limit  cycle  existence,  since  every 
(simple)  friction  limit  cycle  for  this  system  model  meets 
these  conditions.   The  issues  of  sufficient  conditions  and 
minimal  sets  of  necessary  and  sufficient  conditions  are 
discussed.   Finally,  some  examples  applying  these  conditions 
are  presented. 

Necessary  Conditions  for  a  Simple  Limit  Cycle 

In  order  to  simplify  the  analysis,  the  assumption  is 
made  that  the  limit  cycle  which  the  system  exhibits  is 
simple,  in  the  sense  that  it  traverses  each  of  the  four 
state-space  regions  exactly  once  during  each  limit  cycle 
period.   This  is  assumption  (A6)  discussed  in  Chapter  III. 
This  assumption  makes  sense  from  a  physical  viewpoint,  since 
it  is  hard  to  imagine  a  system  where  the  velocity  changes 
sign  more  than  twice  before  returning  to  the  initial  state. 

It  is  clear  that  this  assumption  holds  in  the  2D  case, 
since  the  trajectory  must  go  to  the  left  in  region  I  and  to 
the  right  in  region  II  (see  Figure  1-3  and  1-5)  and  the 
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trajectory  cannot  cross  over  itself.   However,  no  claim  can 
be  made  that  the  following  analysis  is  completely  general, 
due  to  this  assumption.   Therefore,  although  this  is  a 
restriction  of  the  generality  of  the  results,  the  solution 
of  the  equations  requires  the  limit  cycle  meet  the  following 
definition: 

Simple  Limit  Cycle:   A  periodic  orbit  which  traverses  each 
piecewise  linear  region  once  at  most  (at  most  twice  for 
region  III) . 

As  a  consequence  of  this  assumption,  there  will  be 
exactly  four  distinct  switching  instants  during  the  period. 
The  system  will  be  in  a  particular  region,  and,  therefore, 
its  behavior  will  be  governed  by  that  region's  system  model 
equations,  for  T.  seconds,  for  i  =  1,  2,  3,  4.   As  will  be 
seen  below,  some  of  these  switching  periods  may  be  zero. 
Thus  the  total  limit  cycle  period  T  is  the  sum  of  the  four 

T.  's: 

1 

(4.1)      T  =  T   +  T   +  T   +  T 

Note  that  the  subscript  on  the  periods  indicate  order,  and 

not  the  region  being  traversed  during  period  i. 

Since  the  limit  cycle  is  periodic,  any  initial  point 
may  be  chosen  to  start  the  analysis.  Define  x.  to  be  the 
system  state  at  t=0.  Assuming  a  limit  cycle  exists,  this 
initial  state  can  be  chosen  to  be  the  point  in  region  III 
(or  V,  see  Figure  1-3,  p.  6)  where  motion  is  impending.  In 
other  words,  initial  velocity  is  zero,  and  the  initial 
driving  torque  is  negative,  and  sufficient  to  overcome 
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static  friction  and  begin  motion.   These  requirements  on  the 
initial  point  can  be  expressed  as  (recalling  that  y 
represents  the  physical  state  variables,  and  x  the  control 
canonical  states) : 
(4.2a)     y^(0)   =   [  0  ^^   ^^   .  .  .  p^_^       i  ]   x^ 

=   I'  Xq   =   0 
(4.2b)     y^(0)   =   ^'  A  Xq   <   -Lg   (or  =  -Lg  ) 
where  the  notation  j0'  is  used  to  represent  the  n    row  of 
the  x  to  y  transformation  in  (3.8)  (used  to  calculate  the 
velocity  y^  from  the  state  x) ,  the  parameter  Lg  is  the 
maximum  static  friction  before  breakaway,  and  the  matrix  A 
is  the  control  canonical  form  of  the  state  matrix,  in  region 
I.   Note  that  the  A  matrix  from  region  I  (velocity  <  0)  is 
used  since  the  rate  of  change  of  y   is  identically  zero  in 
region  III,  by  design. 

Having  characterized  the  requirements  on  the  initial 
point,  the  limit  cycle  trajectory  can  be  followed  in  region 
I  (which  the  system  transits  after  breaking  free)  using  the 
region  I  model.   The  autonomous  state  variables  can  be 
formed  as  discussed  in  Chapter  III,  by  making  the 
translation: 

(4.3)      ^la   =   ^   -   (-^c/a^l) 


or 


with 


X,   =   X   -   1 
—a     —     — 


1'   =   [  (-Lc/a^^)   0   0 
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where  a  .  is  the  value  in  the  n    row  and  first  column  of  A, 
nl  ' 

and  Lc  >  0  is  the  sliding  coulomb  friction  torque.   This 
translation  is  selected  so  that  the  friction  torque  is 
positive,  and  opposes  the  motion  in  region  I  in  the  negative 
direction.   Note  that  the  parameter  a  ^  cannot  be  zero, 
since  this  would  give  a  zero  column  in  A,  which  signals  the 
presence  of  a  free  integrator. 

The  trajectory  in  region  I  is  then 

(4.4)  x^(t)   =   e^^  x^Q   =   e^^  (x^  -  1) 

so  that  the  trajectory  exits  region  I  at  the  point  where  y 
(velocity)  again  goes  to  zero 

(4.5)  x(T^)   =   x^(T^)   +   1 

AT 

=  e^'^l  (Xq  -  1)   +1 
with  the  requirement  that  y  (T^ )  =  0: 

(4.6)  YnC^i)  =  ^'  ^(^i)  =    ° 

At  this  point,  the  limit  cycle  will  display  two  types 

of  behavior,  depending  on  the  driving  torque  at  t=T- .   If 
the  trajectory  transits  from  region  I  into  region  V,  this 
torque  is  already  sufficient  to  break  loose  the  system  and 
drive  back  in  the  other  direction,  then  T   =  0  and  the 
system  spends  no  time  in  region  V  (in  other  words,  system 
does  not  stick,  but  is  only  momentarily  motionless) .   As 
discussed  in  Chapter  III,  the  region  V  equations  do  not  have 
to  be  applied,  since  the  entry  and  exit  points  for  the 
trajectory  in  this  region  are  the  same.   If  this  also  occurs 
at  the  other  switch  point  (i.e.,  T   =  0),  the  limit  cycle  is 
driven  purely  by  sliding  friction. 
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The  other  case  is  more  complex  to  analyze,  since  the 
trajectory  must  be  calculated  through  all  four  switching 
periods.   If  the  driving  torque  is  insufficient  to  overcome 
the  static  friction  at  t  =  T  ,  the  mechanical  load  stops  for 
a  period  T„  until  the  controller  can  integrate  up  and  break 
it  loose.   The  system  has  come  to  rest  in  region  III  after 
leaving  region  I,  so  system  behavior  is  governed  by  the 
model  for  region  III  during  this  period.   Both  static  and 
sliding  friction  effects  will  affect  this  limit  cycle.   For 
this  case,  the  system  state  evolves  as 

(4.7)  x(t)  =   e^g(^"^^^  x(T^),  T^  <  t  <  (T^  +  T^) 
from  the  initial  point  when  the  velocity  went  to  zero  at 
t  =  T^ .   Note  that  no  translation  is  necessary,  since  the 
friction  term  is  absorbed  into  the  system  model  in  the 
static  condition,  represented  by  state  matrix  A  . 

In  either  the  pure  sliding  friction  case  (T_  =  0)  or 
the  sticky  case,  the  requirements  on  the  system  state  at 
t  =  T-  +  T   is  that  the  driving  torque  is  sufficient  to 
break  loose: 

(4.8)  y^(T^  +  T^)  =  ^'A  x(T^  +  T^)  >  Lg   (or  =  Lg) 
Since  the  sticky  case  involved  the  system  evolving  until 
sufficient  torque  was  developed  to  break  away,  the  system 
will  move  as  soon  as  torque  equals  the  static  friction,  so 

(4.8)  becomes  an  equality;  note  that  the  same  comment 
applies  to  condition  (4.2b). 

There  is  also  the  implicit  requirement  that  the 
velocity  is  zero  when  motion  is  impending  (as  in  condition 
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(4.2)),  but  this  requirement  is  not  placed  explicitly,  since 
it  is  taken  care  of  by  the  model  equations  in  region  III 

(which  force  acceleration  to  be  zero)  plus  the  zero  velocity 
at  t  =  T^. 

The  system  goes  through  the  same  behavior  as  above,  in 
the  opposite  sense,  while  traversing  region  II  and  returning 
to  III  (or  IV) .   The  equations  for  the  trajectory  in  these 
regions,  along  with  the  conditions  on  the  state  at  the 
switching  instants,  are  found  by  first  performing  the 
translation  to  autonomous  coordinates  for  region  II,  as  in 

(4.3): 

(4.9)  x^^  =   X   +   (-Lc/a^^) 
or 

X    =   X   +   1 
—a     —     — 

with 

1-  =   [  (-Lc/a^^)   0   0  .  .  .   0   ] 
as  before.   Note  the  friction  torque  is  negative  in  this 
region  for  L^  >  0 ,  since  the  velocity  is  positive.   The 
trajectory  in  region  II  is  then 

(4.10)  X  (t)  =   e^^^~'^l"'^2^  x^(T,  +  T.) 

a  a   X      z 

=   e^(^"'^l"'^2^  (x(T^  +  T2)  +  1), 
for   T   +  T   <  t  <  T   +  T   +  T 
and  the  trajectory  exits  region  II  at  the  point  where  y 
(velocity)  again  goes  to  zero 

(4.11)  x(T^  +  "^2  ^  "^3^  "   -a^'^1  "^  "^2  "^  "^3^   "   - 

AT 

=  e'^^S  (  x(T^  +  T2)  +   1)   -   1 
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with  the  requirement  that  y  (T   +  T„  +  T  )  =0: 

(4.12)  YnC^i  +  '^2  ^  '^3^  "  ^'  -("^l  "^  "^2  "^  '^3^  "  ° 

At  this  point,  the  system  trajectory  either  immediately 

breaks  free  (trajectory  in  region  IV)  and  leaves  the  region, 
or  it  sticks  (trajectory  in  region  III)  and  follows  the 
trajectory  defined  by 

(4.13)  x(t)  =  e^g^^"^l"^2~'^3^  x(T^  +  ^2  ^  ^3^' 

for   T-  +  T„  +  T_  <  t  <  T-  +  T^  +  T^  +  T^ 
12     3  12     3     4 

from  the  point  when  the  velocity  went  to  zero.   Note  that  no 
translation  is  necessary,  since  the  friction  term  is 
absorbed  into  the  system  model  in  the  static  condition  as  in 
the  case  of  region  III. 

Collecting  the  conditions  (4.5,  4.7,  4.11,  4.13)  for 
the  trajectory  through  all  four  regions,  and  placing  the 
requirement  that  the  end-point  match  the  initial  point,  we 
have  the  following  nonlinear  algebraic  matrix  equation  for 
the  initial  point: 

(4.14)  {I  -  e^s^4   e^'^3   e^s^2   e^"^!}  x^  =  1^ 
with 

1^   =   e^s'^4  (e^'^3  [e^s'^2  {I  -  e^'^l)  1  +  1]  -  1) 

A  T.   AT^   A  T^   AT,  ,  ,   A  T^   AT^   A  T^  , 
=-es4e   3es2e   li+es4e   3es2i 

^   AT.   AT_  ,     AT., 
+  es4e   31-es41 

In  addition,  there  are  four  scalar  nonlinear  equations 

from  the  four  conditions  on  switching  points  (4.2a/b,  4.6, 

4.8,  and  4.12,  but  the  last  one  is  redundant).   Altogether, 

there  are  n+4  (nonlinear)  equations  in  the  n+4  unknowns  x_, 

T  ,  T  ,  T  ,  and  T  .   In  the  pure  sliding  friction  case,  of 
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course,  T   and  T   are  both  zero,  but  conditions  (4.2b)  and 
(4.8)  become  inequalities  and  thus  drop  out,  and  (4.14) 
simplifies  to 

AT     AT 

(4.15)     (I  -  e'^^3    e^'-'l}  x^   =   1^ 
with 

1^   =  e^'^3  [  (I  -  e^'^1}  1  +  1  ]  -  1 
=  -e^'^3  €^"^1  1  +  2  e^'^3  1-1 
so  there  are  n+2  equations  in  the  n+2  unknowns  x  ,  T  ,  and 
T   ,  plus  the  two  inequality  conditions  on  the  switch 
points. 

The  results  presented  so  far  in  this  chapter  can  be 
summarized  as 

Theorem  4-1:   A  simple  limit  cycle  exists  for  the  system 
model  under  study  only  if  the  following  algebraic  conditions 
hold  on  the  variables  x  ,  T  ,  T  ,  T  ,  and  T  : 

Sticking  Limit  Cycle  Case 
(4.14)     {I  -  e^s'^4   e^^3   e^s'^2   e^'^l)  x^  =  1^ 

where      1^  =  e^s'^4  (e^'^3  [e^s'^2  {I  -  e^"^!}  1  +  1]  -  1) 

A  T^   AT^   A  T^   AT^  ,  ^   AT.   AT_  ^A^T_  , 
=-es4e   3es2e   ll+es4e   3es21 

A  T^   AT^  ,     A  T^  , 
+  es4e   3i-es4i 

(4.12)     y^(T^  +  ^2  +  T3) 

=  [  0  ^0   ^1   •  •  •   ^n-3   ^  ^  -^^-^  +  T^  +  T.,) 


(4.2b)     Y^{0)       =      I'  A  Xq  =  -Ls 

(4.6)      Yn^^i)  =  ^'  ^(T^)  =  0 

(4.8)      yj^(T^  +  T^)  =  £'A  x(T^  +  T2)  =  Lg 
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Pure  Sliding  Limit  Cycle  Case  (T2  =  T3  =  0) 


(4.15)     {I  -  e^'^3  e-^^l)  x^   =   1^ 

AT  AT 

Where     1^   =  e   3  [  {I  -  el)  1  +  1  ]  -  1 

=  _eAT3  ^AT^  ^  ^  ^  e^^S  1  -  1 

(4.2a)     y^(0)   =   [  0  ^^   ^^   ...   ^^,3   i  ]   x^ 

=  ^'  Xq   =   0 

(4.6)      Yn^^l^  =  ^'  ^(^1^  =  ° 

In  other  words,  the  conditions  above  are  necessary  for  the 

existence  of  a  simple  limit  cycle.   Condition  (4.12)  was 

included  in  the  sticking  case  instead  of  the  equivalent 

condition  (4.2a)  to  obtain  a  more  symmetric  form  of  the 

conditions. 

Note  that  the  conditions  requiring  velocity  =  0  could 
be  applied  at  either  t  =  T^  or  t  =  T^  +  T2  (equation  (4.6)), 
and  at  either  t  =  0  or  t  =  T^  +  T2  +  T3  (equation  (4.2a)). 
This  is  a  consequence  of  the  fact  that  velocity  remains  zero 
between  these  times. 

Proof;   The  derivation  presented  above  is  a  step-by-step 
solution,  through  each  region,  for  the  limit  cycle 
trajectory  that  was  assumed  to  exist.   Since  the  conditions 
stated  in  the  theorem  were  derived  directly  from  the 
trajectory  solutions  as  seen  above,  and  since  a  simple  limit 
cycle  by  definition  must  traverse  those  regions  as  assumed, 
the  conditions  are  necessary,  and  the  proof  is  complete. 

Note  that  a  non-simple  limit  cycle,  should  any  exist, 
would  not  have  to  meet  the  necessary  conditions.   Thus,  if  a 
solution  to  the  necessary  conditions  exists  for  a  system,  a 
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limit  cycle  may  exist;  if  no  solution  exists,  no  simple 
limit  cycle  exists,  but  a  non-simple  one  is  not  ruled  out! 
A  complete  theory  requires  an  examination  of  the  possibility 
of  more  complex  limit  cycles. 

It  is  interesting  to  note  that  the  same  development 
could  be  performed  using  physical  state  variables.   This 
would  result  in  replacing  the  control  canonical  A  matrix  in 
the  equations  by  the  original  system  matrix,  eliminating  the 
need  for  the  vector  §_,    and  replacing  1  by  the  translation  in 
equation  (3.12).   The  form  used  here  clarifies  the  role  of 
system  poles  and  zeroes,  however,  as  seen  by  the  following 
discussion. 

The  definition  of  the  initial  state  Xq,    equation 
(4.14),  is  completely  determined  by  the  poles  (for  fixed 
switching  periods) ,  since  the  A  matrix  contains  the 
characteristic  polynomial  information  only.   Now,  observing 
that  the  transfer  function  from  the  friction  input  to  the 
velocity  as  an  output  is 

yn(s)  /  l(s)  =  [0  0  .  .  .  0  1]  Ti2  (sI-A)-l  b 
(this  holds  since  y  (s)=T-l2X(s)  ,  and  x(s)  =  (sI-A)  "^bl  (s)  ), 
with  b'  =  [0  0  .  .  .  1]  ,  it  is  clear  that  the  transfer 
function  zeroes  are  defined  by  the  bottom  row  of  Ti2'  that 
IS,  the  vector  §_.      Thus,  varying  the  zeroes  only  would 
change  the  ^  vector  only  (in  such  a  way  that  poles  of  the 
closed  loop  system  are  unaffected),  Xq  would  be  unchanged, 
and  the  four  (or  two)  switching  conditions  alone  would  be 
affected.   The  control  form  representation  clarifies  this 
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relationship,  which  suggests  an  interesting  line  for  further 
research. 

Illustrative  Examples 
Example  IV-1;   Two-Dimensional  (2D)  System  with  Sliding 

Friction  Limit  Cycle 

This  is  the  same  as  Example  I-l  of  Chapter  I;  the 
piecewise  linear  method  is  the  same  as  standard  phase  plane 
methods  for  this  2D  case.   The  state  variable  model  of  the 
system  is 
(4.16)     dyi/dt   =   y2 

dy2/dt   =   -(K/J)  yi   -   (B/J)  Y2      -    Lf/J 
where  Lf  is  the  friction  torque,  y^  is  position,  and  Y2    is 
velocity.   The  analysis  performed  below  assumes  K,  J  >  0, 
although  B  can  be  negative  (the  other  cases  can  be  solved  by 
similar  methods,  but  are  left  out  for  brevity) .   Chapter  I 
presents  more  information  about  the  system,  including 
results  of  a  previous  solution  by  phase  plane  methods. 

A  second-order  system  cannot  possibly  have  a  limit 
cycle  with  sticking,  since  once  the  system  sticks,  it  is  at 
an  equilibrium  point  (there  is  insufficient  torque  to  break 
free,  and  no  controller  integrator  to  ramp  up) .   Therefore, 
we  need  only  examine  the  necessary  conditions  for  the 
sliding  case.   In  addition,  it  will  be  assumed  for  this 
example  that  the  limit  cycle  is  symmetric  (this  property  is 
proved  for  2D  systems  in  Chapter  V) .   Applying  the  equations 
developed  in  this  chapter,  the  initial  condition  can  be 


60 


found  from  equations  (4.5)  (which  (4.15)  simplifies  to  using 
the  symmetry  assumption)  and  (4.6): 

AT 

(4.17)  x(T^)  =  e^'^l  (Xq  -  1)   +1   =   -xo 

(4.18)  Ys^^i)  =  ^'  ^C^i)  =  ^2(Ti)  =  0 

where  use  is  made  of  the  fact  that  ^'=[01]  (physical 
state  model  already  in  canonical  form),  and  x(T^)  =  -Xq  (by 
symmetry) . 

Making  the  assumption  that  K  is  large  enough  (or  B  is 
small  enough)  that  the  system  poles  are  complex,  the  matrix 
exponential  for  this  2D  system  can  be  evaluated  as 


(4.19) 


At  -at 

e        =   e 


cos)3t+(a/j0)sin/3t       (l//?)sin;3t 
-(K/J/3)sir\l3t      cos)3t-(a//3)sin/3t 


where      a   =   B/(2J),  0   =    (1/2)  (  (4K/J)  -  (B/J)2  }  ^^ . 
Letting  Xq  =  [a  0]'  (since  X2=0) ,  and  noting  1  =  [Lq/K  0]', 

equations  (4.17)  and  (4.18)  become: 

— aT 
(4.20)     X2(Ti)  =  0  =  -[a  -  (Lq/K)  ]  e    1  (K/J/3)  sin/3Ti 

which    implies   T^    =   7r//3    (note   a    >    (Lq/K)    since    otherwise   Xg 

is  an  equilibrium  point) . 

Theoretically,  T^   can  also  be  any  integer  multiple  of 

-at    . 
7r//3.       However,    noting   that   X2  (t)    =   velocity   =   C   e         sin/3t, 

any  multiple  greater  than  one  would  entail  sign  changes  in 

the  velocity  for  t<T2  (i.e.,  trajectory  leaves  region  II). 

This  is  an  example  of  a  solution  to  the  necessary  conditions 

that  is  not  a  valid  limit  cycle;  see  discussion  on 

sufficient  conditions  below. 
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The  other  component  of  xq  can  be  evaluated  from  (4.17) 
as 

(4.21)  Xi(Ti)  =  -a 

=  [a  -  (Lc/K)  ]  e"'^'^l(cos)3Ti+(a/^)sin/3Ti)  +  (Lq/K) 
=  [a  -  (Lc/K)  ]  e""^/^  costt  +  (Lq/K) 
Therefore 

(4.22)  a  +  (Lc/K)  =  [a  -  (Lc/K)]  e~^'^^^ 
The  initial  condition  can  now  be  evaluated  as 

(4.23)  a  =  (Lc/K)  {  (m+1)  /  (M-1)  )  ,  M  =  e""""/^ 

which  is  the  result  presented  in  Chapter  I.   Since  Xq  is  in 
region  IV,  a  must  be  greater  than  zero.   This  implies  ii>l , 
requiring  a<0,  or  B<0. 

As  the  solution  from  Chapter  I  predicted,  a  sliding 
type  of  limit  cycle  exists  for  every  case  with  negative 
damping  (B<0)  (if  the  position  feedback  gain  K  is 
sufficiently  large  to  give  the  linear  system  a  complex  pair 
of  poles) .   The  initial  conditions  can  be  expressed  in  terms 
of  the  system  parameters  by  the  equations  from  Chapter  I: 

(1.2)  Xio  =  (Lc/K)  [(M+1)/(M-1)],    X20  =  0 
where 

(1.3)  M  =  exp[-B7r/(2J)3)  ],  p    =    h    [4K/J  -  b2/j2]^ 
Refer  to  Figure  1-4  for  a  plot  of  amplitude  (which  equals 
Xjo)  3S  a  function  of  gain  K  for  the  case  J=Lc=l,  B=-l. 
Example  IV-2;   Three-Dimensional  (3D)  System  with  Both 

Sliding  and  Sticking  Limit  Cycles 
The  second  example  system  to  be  considered  (same  as 
Example  1-2)  is  defined  by  the  differential  equations 
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(1.4) 
where 
(1.5) 


dy/dt   =   Ay   +   b  L. 


0 

0 

-K, 


1 

0 

-K, 


0 

1 
-B 


1      "2 
with  b'  =  [0  0  -1],  L^  =  friction  torque,  and  where  the 

state  variables  are  y^  =  position,  y  =  velocity,  and 

y^  =  compensator  integrator.   Figure  1-6  shows  a  block 

diagram  of  this  system.   The  results  presented  in  Chapter  I 

are  now  to  be  derived,  based  on  the  equations  developed  in 

this  chapter. 


Case  I;   One  Stable  Pole  and  One  Imaginary  Pole  Pair  (Only 
Sticking  Limit  Cycle  Exists) 

As  in  the  example  in  Chapter  I,  let  us  first  set 
B  =  K^  =  K^  =  1,  and  set  the  sliding  friction  torque  to  1.0 
and  sticky  friction  (breakaway  torque)  to  1.2. 

The  equations  (4.2a/b,  4.6,  4.8)  defining  the  switching 
periods  and  initial  state  (4.14)  can  be  set  up  and  solved 
for  this  specific  case,  in  order  to  demonstrate  the 
analytical  calculation  of  the  limit  cycle  trajectory.   The 
eigenvalues  of  this  system  are  -1  and  ±jl.   The  matrix 
exponential  of  A  is  required  for  the  equations,  and  can  be 
found  (by,  for  example,  diagonalizing  A)  as 


(4.24) 


2e 


AT 


-T   .  .  -T 

e   +sinT+cosT  2sinT  e   +sinT-cosT 

-T   .  -T 

-e   -sinT+cosT  2cosT  -e   +sinT+cosT 

-T   .  .  -T 

e   -sinT-cosT  -2sinT  e   -sinT+cosT 
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The  equations  to  be  applied  are 
(4.25)     --     -'^'^- 
(4.26) 
where  a  symmetric  limit  cycle  was  assumed  (breakaway  point 


x^  =  e   1  (Xq  -  i)  +  i 

A  T^ 
x^  =  e  s  2  x^  =  -Xq 


at  t  =  T^  +   T-    is  at  -Xf^)  ,  and  therefore  the  simpler 
equations  were  used  in  place  of  (4.14).   Note  that  A   is  the 
A  matrix  with  the  bottom  row  set  to  zero,  hence  the  matrix 
exponential  in  the  sticking  region  is 


(4.27) 


e^sT2 


1    T2   T2V2 

0    1     T2 

0    0      1 

In  addition,  the  breakaway  condition  that  the  acceleration 
torque  equal  the  static  friction  is  required: 


(4.28) 


^'Ax 


0 


dX3(0)/dt  =  -Lg 


Noting  that  the  translation  vector  1'  =  [  10  0  ],  letting 
X-  =  [  a  b  0  ]  (since  x  .  must  be  zero) ,  and  letting 


X]^  =  [  c  d  0  ]  ,  we  obtain 


-T 


(4.29)  c  -  1  =  (1/2) (a-1) (e    +  sinT  +  cosT)  +  bsinT 

-T 
d  =  (1/2)  (a-1)  (-e    -  smT  +  cosT)  +  bcosT 

0  =  (1/2)  (a-1)  (e~   -  sinT  -  cosT)  -bsinT 

from  application  of  equation  (4.25)  and  the  definition  of 

AT 
the  matrix  exponential  e   , 

(4.30)  c  +  d  T2   =   -a 
d   =   -b 

from  application  of  equation  (4.26)  and  the  definition  of 
e  s  2,  and 


(4.31) 


a  +  b   = 
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from  application  of  equation  (4.28).   By  eliminating  c  and  d 
from  these  equations  and  performing  additional  algebra,  they 
simplify  to 
(4.32)     b  =  Lg  -  a 

T^  =  [(a-l)e"'^l  +  (a+i)]  /  (Lg  -  a) 

-T 
a  [e   1  +  sinT   -  cosT  ] 

=  e~'^l  +  (21^^  -  l)sinT^  -  COST 

-T 
a  [-e   1  -  sinT   -  cosT   -2] 

-T 
=  -e   1  -  (2  Lg  -  l)cosT^  -  sinT^  -  2L 

The  first  two  equations  define  b  and  T   in  terms  of  a  and 

T^,  while  b  and  T2  were  eliminated  from  the  last  two,  which 

can  be  used  to  simultaneously  solve  for  a  and  T  . 

Eliminating  a  and  after  some  algebra,  the  last  two  equations 

form  a  nonlinear  algebraic  equation,  whose  zeroes  are 

potential  limit  cycle  solutions: 

(4.33)     f(T^)  =  (e~^l  -  1) (1+cosT^)  -  (e"^l  +  l)sinT^  =  0 

The  plot  of  this  function  in  Figure  (IV-1)  and  analysis 

shows  that  there  are  zeroes  at  odd  multiples  of  ir ,    and  also 

near  37r/2  (+  2n7r,  n  =  0,  1,  ...).   The  solution  at  n   results 

in  a  valid  limit  cycle  solution  where  x  •  =  [  1   0.2   0  ], 

T^  =  n ,    and  T^  =  10,  yielding  an  overall  limit  cycle  period 

of  approx.  26.3  seconds,  matching  the  simulation  results. 

The  second  potential  solution  is  at  T   =  4.73  (approx.  37r/2, 

obtained  by  numerical  solution  of  f(Ti)=0),  and  is  invalid, 

since  it  results  in  a  negative  value  for  T   (of  -12.214). 

Additional  solutions  are  also  invalid,  for  the  same 

reason  as  in  Example  IV-1.   Examination  of  the  matrix 
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exponential  and  the  resulting  formula  for  the  velocity  state 
shows  that  the  velocity  changes  sign  approximately  every  n 
seconds.   Therefore,  the  trajectory  using  these  longer 
values  for  T^  leaves  region  II  before  the  switching  time,  so 
it  is  invalid  (see  next  section) . 

Therefore,  the  solution  of  the  piecewise-linear 
equations  results  in  a  valid  limit  cycle  solution  that 
matches  simulation  results.   Chapter  I  contains  plots  of  the 
limit  cycle  trajectory. 

Case  II:   One  Stable  Pole  and  One  Unstable  Pole  Pair  (Both 
Sticking  and  Sliding  Limit  Cycles  Exist) 

If  the  same  example  system  is  used  with  B  =  0.9,  and 
K^  =  K^  =  1  still  (case  II  in  Chapter  I) ,  the  closed-loop 
eigenvalues  (of  the  linear  portion  of  the  system)  move  into 
the  right-half  plane,  to  .026  ±  j  1.024,  while  the  other 
pole  is  at  -.9524.   As  stated  in  Chapter  I,  there  is  still  a 
sticking  limit  cycle  solution  close  to  that  of  Case  I 
(Figures  1-7  through  1-9,  confirmed  by  simulation).   A 
similar  analysis  to  that  presented  above  can  be  used  for 
this  case,  so  the  details  are  omitted. 

Note,  however,  that  numerical  methods  may  be  required 
to  evaluate  the  matrix  exponential  (due  to  the  tediousness 
of  the  calculations;  case  I  was  a  particularly  simple  form) ; 
this  forces  the  use  of  iterative  methods  to  calculate  the 
limit  cycle  parameters  (see  example  of  this  method  below) . 

The  second  limit  cycle  solution,  as  discussed  in 
Chapter  I,  is  of  the  sliding  type.   Simulation  having  shown 
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one  with  a  half-period  T   approximately  equal  to  3 . 1 
seconds,  an  iterative  numerical  solution  of  the  equations 
for  the  symmetric,  sliding  case  was  performed,  as  follows. 

The  equation  for  this  case  is  obtained  from  (4.5), 
where  again  the  assumption  is  made  that  the  limit  cycle  is 
symmetric: 

AT 

(4.34)  X-L  =  e^^^l  (x^  -  1)  +  1 

which  can  be  solved  as 

(4.35)  Xq   =  {I  +  e^'^1}-!  {I  -  e^^l)  1 

A  computer  program  calculated  values  for  the  initial 
condition  Xq,  given  trial  values  of  Tj,  and  the  process  was 
repeated  until  the  velocity  initial  condition  Xq^"^^  ^  °- 
This  was  not  difficult,  since  an  approximate  starting 
condition  was  available  from  the  simulation  results.   Given 
a  T]^  that  resulted  in  a  zero  initial  velocity,  the  entire 
initial  condition  could  be  defined  from  (4.35). 

This  yielded  a  solution  at  T   =  3.1455  seconds,  and 
Xq  =  i^O  ^  [-0.283   12.4   0.].   The  solution  is  valid  since 
the  breakaway  condition  is  also  satisfied  (torque  at  zero  is 
greater  than  LS) .   Figures  I-IO  and  I-ll  show  some  views  of 
the  limit  cycle  for  this  case.   An  appendix  contains 
computer  code  that  was  used  to  perform  this  iterative 
calculation.   The  code  applies  to  the  sliding  case  only,  and 
would  have  to  be  modified  for  the  general  case. 
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Exact  (Necessary  and  Sufficient)  Conditions 
Theorem  4-1  presented  some  conditions  that  were 
necessary  for  the  existence  of  a  simple  limit  cycle.   That 
they  are  not  sufficient  is  seen  by  the  results  in  the 
Examples,  where  solutions  were  found  that  met  all  of  the 
necessary  conditions,  yet  were  not  valid  limit  cycles  (the 
switching  period  T2  was  less  than  zero,  or  velocity  changed 
sign  during  a  switching  period) . 

In  order  to  expand  the  previous  conditions  into  a  set 
that  is  also  sufficient,  the  concept  of  a  consistent 
solution  is  useful: 

Definition;   A  consistent  solution  is  a  solution  of  the 
equations  in  Theorem  4-1  that  meets  the  assumptions  about 
the  region  containing  the  trajectory  during  each  switching 
period.   These  assumptions  are 
(Sticking  Case) 

(1)  Trajectory  in  region  I,  0  <  t  <  T^ 

(2)  Trajectory  in  region  III,  T]^  <  t  <  T^  +  T2 

(3)  Trajectory  in  region  II, 

(4)  Trajectory  in  region  III, 

T;]L  +  T2  +  T3  <  t  <  Tj  +  T2  +  T3  +  T4 
(Pure  Sliding  Case) 

(1)  Trajectory  in  region  I,  0  <  t  <  T^ 

(2)  Trajectory  in  region  V,  t  =  T^ 

(3)  Trajectory  in  region  II,  T^  <  t  <  T^  +  T3 

(4)  Trajectory  in  region  IV,  t  =  T^  +  T3 
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A  check  of  the  consistency  of  a  solution  is 
straightforward,  based  on  the  definitions  of  the  various 
regions  presented  in  previous  chapters:   the  sign  of  the 
velocity  is  examined  for  assumptions  (1)  and  (3),  and  the 
magnitude  of  the  torque  for  (2)  and  (4) .   Given  this 
definition,  a  set  of  necessary  and  sufficient  existence 
conditions  can  now  be  presented: 

Theorem  4-2;   A  simple  limit  cycle  exists  in  the  system 
under  study  if  and  only  if  a  consistent  solution  of  the 
Theorem  4-1  equations  exists.   In  addition,  this  limit  cycle 
has  the  exact  properties  (switching  periods,  initial 
condition,  and  trajectory)  defined  by  the  corresponding 
solution  to  the  equations. 

Proof:   (Necessity)   Suppose  a  simple  limit  cycle  exists 
with  the  given  properties.   Theorem  4-1  (more  exactly,  the 
derivation  in  the  first  part  of  this  chapter)  demonstrated 
by  analysis  of  this  limit  cycle  that  the  listed  equations  in 
fact  apply  to  the  trajectory. 

In  addition,  this  derivation,  along  with  the  assumption 
that  the  limit  cycle  is  simple,  shows  that  the  solution  must 
meet  the  consistency  conditions,  i.e.,  it  must  be  in  the 
appropriate  region  of  the  state  space  during  each  switching 
period.   Therefore  this  condition  is  also  necessary. 

(Sufficiency)   Now  suppose  a  consistent  solution  of  the 
equations  exists.   A  limit  cycle  can  be  generated,  with  the 
same  properties  as  defined  in  the  solution,  by  starting  at 
Xq  and  following  around  to  the  initial  point  again.   At  each 
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step  of  the  process,  the  consistency  conditions  require  the 
trajectory  to  be  within  a  given  region,  so  that  the  system 
equations  for  that  region  apply.   The  equations  listed  in 
Theorem  4-1,  derived  from  the  application  of  those  system 
equations,  then  show  that  the  trajectory  reaches  the  next 
switching  point  at  the  appropriate  time  and  state. 
Therefore,  a  limit  cycle  must  exist  with  those  exact 
properties,  and  the  proof  is  complete. 

Although  this  theorem  defines  a  set  of  exact 
conditions,  i.e.,  they  are  equivalent  (necessary  and 
sufficient)  to  the  existence  of  the  limit  cycle,  it  is 
appropriate  to  ask  if  these  conditions  are  minimal.   In 
other  words,  are  there  redundancies  in  the  set  of 
conditions?   Is  there  a  simpler  set  of  conditions  that  are 
still  exact?   For  example,  it  may  be  sufficient  to  merely 
require  switch  periods  greater  than  zero,  thus  eliminating 
cases  such  as  found  in  Example  IV-2 .   These  questions  are 
open,  and  may  be  subjects  of  further  research. 


CHAPTER  V 
LIMIT  CYCLE  SYMMETRY  AND  STABILITY 

This  chapter  uses  the  results  of  chapter  IV  and  other 
known  facts  about  the  friction  nonlinear ity  to  examine  the 
behavior  and  characteristics  of  friction  limit  cycles.   The 
model  of  the  system  developed  previously  is  used  to  analyze 
the  stability  of  limit  cycles  predicted  by  the  equations  in 
chapter  IV.   In  addition,  the  relationships  between  the 
oddness  of  the  system  differential  equations,  symmetry,  and 
uniqueness  of  limit  cycles  is  explored. 

Stability  of  Predicted  Limit  Cycles 

Once  a  limit  cycle  has  been  predicted  by  the  solution 
of  the  equations  derived  in  Chapter  IV,  it  is  natural  to  ask 
about  its  relationship  to  the  global  phase  portrait.   Since 
the  piecewise  linear  model  used  here  is  exact,  the  exact 
limit  cycle  period,  amplitude,  and  trajectory  can  be  found 
immediately  once  an  initial  point  on  the  limit  cycle  is 
known.   It  is  more  difficult  to  determine  stability  and 
other  information  about  the  phase  portrait,  however,  since 
it  can  require  extensive  simulation. 

However,  local  stability  of  the  limit  cycle  can  be 
determined  by  the  standard  method  of  linearization  of  the 
periodic  flow  map  at  the  initial  point  previously  found.   In 
other  words,  an  approximate  model  of  the  trajectories  near 
the  limit  cycle  initial  point  can  be  determined  by 
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linearization.   This  analysis  provides  information  about  the 
asymptotic  stability  of  a  closed  orbit  solution  as  discussed 
in  Hirsch  and  Smale  (1974),  chapter  13  (if  more  information 
on  orbital  stability  analysis  is  of  interest,  this  reference 
provides  a  good  description,  and  was  quite  useful  in  the 
development  of  the  results  in  this  chapter) . 

Note  that  the  model  provided  by  this  analysis  is  not  a 
continuous-time  representation  of  trajectories  near  the 
cycle;  instead,  this  discrete-time  model  shows  the  deviation 
from  the  initial  point  and  how  it  evolves  after  each  cycle 
period.   Specifically,  an  initial  sufficiently  small 
deviation  from  the  initial  point  x   is  assumed  (so  the 
trajectory  starts  at  Kq+SXq    )  and  the  trajectory  traced  once 
around  the  limit  cycle.   The  new  deviation  from  the  initial 
point  when  the  trajectory  completes  the  cycle  is  related  to 
the  initial  deviation  by  the  equation 
(5.1)      6x^  =   Y  5Xq   +   higher  order  terms 
where  6x^   is  the  new  deviation  from  the  initial  point 
(Kq+Sx^    is  the  point  on  the  trajectory  after  one  cycle) ,  and 
Y  is  the  transition  matrix,  representing  the  first-order 
evolution  of  the  deviations  from  the  initial  point  over 
time.   Eigenvalues  of  the  transition  matrix  determine 
stability,  while  eigenvectors  determine  the  phase  portrait 
near  the  limit  cycle. 

This  analysis  is  attempted  below;  however,  because  of 
the  discontinuities  in  the  system  nonlinearity  the  analysis 
is  difficult.   Therefore,  the  analysis  is  accomplished,  as 
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in  the  derivation  of  Chapter  IV,  by  solving  one  region  at  a 
time,  and  then  pasting  the  trajectories  together.   The  first 
case  analyzed  is  that  where  the  predicted  limit  cycle  has  no 
sticking  (sliding  friction  only) .   The  case  with  sticking 
during  the  cycle  is  then  analyzed  by  building  on  the  first 
case. 
Case  1:  Limit  Cycle  with  no  Sticking 

Assume  an  initial  point  x.  that  meets  the  equations  in 
Chapter  IV  with  no  sticking.   In  this  case,  the  first 
switching  period  is  equal  to  the  constant  T^  in  chapter  IV, 
the  second  switching  period  is  T  ,  and  the  constants  T   and 
T  ,  representing  the  periods  of  stickiness,  are  zero. 

A  deviation  5x^  is  assumed  from  the  initial  point  of 
the  limit  cycle.   We  wish  to  limit  deviations,  however,  so 
that  the  perturbed  initial  point  is  still  contained  in 
region  IV.   This  is  necessary  in  order  to  apply  the  Chapter 
IV  equations,  since  the  initial  point  to  which  the  equations 
apply  is  assumed  to  be  in  region  IV.   Limiting  the  deviation 
in  this  fashion  results  in  no  loss  of  information  about  the 
limit  cycle  stability,  however;  (n-1)  of  the  eigenvalues  of 
the  periodic  flow  map  can  still  be  determined,  while  the  n 
eigenvalue  (that  applies  to  deviations  out  of  region  IV,  and 
along  the  limit  cycle)  would  be  zero,  as  will  be  seen  in  the 
examples  below.   Note  that  the  derivative  of  the  flow  after 
one  orbital  period  has  an  eigenvalue  of  unity  (Hirsch  and 
Smale,  1974,  p.  277),  but  the  Poincare  map  has  a  zero 
eigenvalue  (Hirsch  and  Smale,  chapter  13,  section  3). 
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In  order  to  limit  the  deviation  Sx      such  that  the 
initial  point  remains  in  region  IV,  <5x   can  be  constructed 
as 

(5.2)      5Xq  =  T^^-l  i^'  6yQ^ 

where  T^^  is  the  composite  linear  transformation  from 
Chapter  III  used  to  convert  between  canonical  and  physical 
state  variables  (y  =  T^^  ^)  '  ^j-   is  an  n  x  (n-1)  matrix  (I 
is  (n-1)  X  n)  used  to  force  the  n^^  element  of  the  deviation 


vector  5y  to  be  zero: 


(5.3) 


I   = 

r 


10  0. 
0  10. 
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and,  finally,  (Sy^^  is  the  (n-1) -vector  of  deviations  in 
physical  state  variables.   This  construction  forces  the  <Sx 
vector  to  remain  in  region  IV,  since  the  change  in  initial 
velocity  {6y_q(^)  ,    the  n^^  element  of  the  physical  state 
vector)  is  zero.   Although  the  derivation  below  will  be  done 
in  canonical  variables,  this  construction  (5.2)  will  be  used 
at  the  end  to  complete  the  calculation  for  physical 
variables,  and  force  the  deviation  to  have  zero  initial 
velocity. 

Starting  at  the  initial  point  Xq+<5x  ,  the  trajectory  is 
determined  by  the  piecewise  linear  methods  as  in  chapter  IV. 
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The  trajectory  in  region  I  is  found  using  autonomous  state 
variables: 

(5.4)  x^(t)  =  x(t)  -  1 
so  the  trajectory  is 

At  At 

(5.5)  ^^(t)  =  e"^""  x^(0)  =  e''''  [Xq  +  Sx^    -    1] 

At 

x(t)  =  e""^  [Xq  +  5Xq  -  1]  +  1 
For  a  sufficiently  small  deviation  <Sx  ,  the  trajectory 
reaches  region  V  (velocity  goes  to  zero  preparatory  to 
reversing  direction) ,  but  not  necessarily  in  exactly  T 
seconds.   Denoting  the  change  in  this  period  is  ST,    the 
velocity  goes  to  zero  when  the  trajectory  reaches 

(5.6)  x(T^  +  5T)  =  e^^'^l"^'^'^^  [x^  +  .Sx^  -  1]  +  i 

The  quantity  of  interest  is  the  deviation  of  this  point  from 
the  point  where  the  original  limit  cycle  reaches  region  V 

AT 

(5.7)  x^  =  e^-^l  [Xq  -  1]  +  1 

(5.8)  6x^  =  X(T^+5T)  -  x^ 

,  A6T    ^,   AT^  ,      ,,  ^   A(T^+(5T)  ^ 
=  (e     -  1}  e   1  [Xq  -  1]  +  e  ^  1      Sx^ 

Note  that  the  first  term  shows  the  dependence  on  61,    since 

it  is  zero  if  6T  is  zero,  while  the  second  term  shows  the 

dependence  on  <5x  . 

Now  that  the  exact  value  for  5x^  is  determined,  this 

nonlinear  formula  must  be  linearized  to  determine  local 

behavior.   One  factor  that  causes  difficulties  is  that  6T 

depends  on  <5x^,  which  complicates  the  calculation  for 

(5.9)  Y  -  (d(5x^)/d(6xQ)}|^^^Q 

where  Y  is  the  desired  transition  matrix  defined  above,  and 
the  derivative  is  evaluated  for  6x     =   0.   The  dependence  of 
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5x^  on  both  ST   and  5x   explicitly  is  handled  using  the  chain 
rule. 

Note  that  since  x   is  not  a  function  of  <5x   or  <5T, 

(5.10)  [d/d(5xQ)]  {Sx^}    =       [d/d(5xQ)]  {x(T^+5T)  -x^) 
=  [d/d(5xQ)]  {x(T^+5T)} 

=  [d/d(5xQ)]  {e^^'^l^'^'^)  [Xq  +  5Xq  -  1]} 
where  the  term  in  x(T  +6T)  equal  to  1  was  also  dropped  since 
it  does  not  depend  on  «5x  .   Using  the  chain  rule,  this 
derivative  is 

(5.11)  d(5x^)/d(<5xQ)  =  D(6x^)/D(<SXq) 

+  [D(<5x^)/D(5T)]  [d(6T)/d(6xQ)  ] 
where  the  symbol  "D"  indicates  a  partial  derivative.   The 
first  term  is 

(5.12)  D(5x^)/D(5Xq)  =  e^^'^l^'^'^^ 

=  e^^l 
when  evaluated  at  <5x  =0  (hence  ST=0)  .   The  partial 
derivative  in  the  second  term  can  be  evaluated  as  follows. 
First  define 

AT 

m  =  e  ^1  [Xq  +  5Xq  -  1] 
Then,  by  using  (5.10),  the  desired  partial  derivative  is 

(5.13)  D(5x^)/D(<5T)   =  [D/D(5T)]  (e^'^'^m) 

-   A(5T 
=  A  e     m 

=  A  e'^^l  [Xq  -  1] 

when  the  equation  for  m  is  substituted  back  in,  and  the 

derivative  is  evaluated  at  <5x  =0  (6T=0).   The 

differentiation  step  can  be  verified  by  expanding  the 
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exponential  into  a  power  series,  and  evaluating  term  by 
term. 

The  evaluation  of  the  [d(5T)/d(5x  ) ]  factor  in  the 
second  term  requires  the  definition  of  a  functional 
relationship  between  ST   and  Sx    .      We  know  from  the 
derivation  of  x(T   +  5T)  above  that  the  velocity  goes  to 
zero  at  this  point  (trajectory  reaches  region  V) .   As  in  the 
analysis  in  Chapter  IV,  this  can  be  stated  as 

(5.14)  §_'    x(T^  +  5T)  =  0 

=  £•  {e-^^^l^^^^  [Xq  +  5Xq  -  1]  +  1} 
=^.  (e^(Ti+^T)  [x^  .  S.^    -1]} 
since  ^'1  =  0.   This  formula  implicitly  defines  the 
relationship  between  ST   and  Sx    ,    in  that  we  have  a  function 
of  the  form  f(5T,5x  )  =  0.   The  implicit  function  theorem 
(Rudin  (1976),  pp.  223-8)  can  be  used  to  evaluate  the 
desired  derivative  as 

(5.15)  d(<ST)/d(5xQ)  =   -  [Df/D(5T)]"^  [Df/D(5xQ)] 
Evaluating  these  factors, 

(5.16)  Df/D(<5T)  =  [D/DiST)]     {§.'  {e^'^'^l'^^'^^     [x^+Sx^-l])  ) 

=    [D/D(5T)]  (£'  e^"^^  m} 

oi  A   A6T 
=  ^'  A  e    m 

with  m  defined  as  above  (to  see  this,  again  expand 

exponential  into  series  and  evaluate  term  by  term) .   The 

second  factor  is 

(5.17)  Df/D(5Xo)    =    [D/D(5Xq)]     {£ '  (e^  ("^l^*^^^     [Xq+^Xq-I  ]  )  ) 
=    [D/D(5Xq)]     {^'    e^^V"^"^^     5Xq) 

=  ^.    e^^^l+^T) 
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The  derivative  from  the  implicit  function  theorem,  (5.15), 
can  now  be  evaluated  as 

AT  —1      AT 

(5.18)  d(5T)/d(5XQ)  =  -{^'A  e^'^l  [Xq-I^  >   [^'^       1) 

where  the  factors  are  again  evaluated  at  i5x_=0  (<5T=0)  . 
Finally,  combining  (5.12,13,18)  into  (5.11),  we  get  the 
total  derivative  representing  the  linearized  periodic  map 

(5.19)  Y  =  {d(5x^)/d(5xQ) }  I^^^Q 

=  e^^l  -  {A  e^^l  [Xq-1]}  (^'A  e^'^1  [x^-l]  }"^(£' e^^l } 
Note  that  the  factor  to  be  inverted  is  a  scalar.   This 
expression  can  be  simplified  by  noting  that 

AT 

(5.20)  x^  =  el  [Xq  -  1]  +  1 
so  substituting  into  (5.19)  gives 

(5.21)  Y  =  e'^'^1  -  {A  [x^-1]}  {^'A  [x^-l  ]  }  "^  {£' e^^l } 

AT  AT 

=  e^'-'l  -  {A  [x^-1]  I'e^'-'l)  /  {£'A  [x^-1]  } 
The  linearized  flow  map  for  the  entire  cycle  must  now 
be  constructed,  using  this  equation  (5.19)  for  the 
linearized  flow  map  for  a  half-cycle.   Note,  first  of  all, 
that  the  linear  transformation  representing  the  linearized 
flow  map  for  the  cycle  is  the  product  of  the  linear 
transformations  of  the  two  halves.   This  is  a  consequence  of 
the  chain  rule,  since 

(5.22)  <Sx^  =  fi(5x^)  =  f;^(f2(<5xQ)) 
d(5x^)/d(5xQ)  =  {d(<5x^)/d(5x^)  )  (d{Sx^) /d{Sx^)  } 

The  second  of  these  two  linear  transformations  has  already 
been  derived;  the  first  must  be  found  and  the  linear  map  for 
the  entire  cycle  found  by  multiplication. 
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Assume  first  of  all  that  the  limit  cycle  is  symmetric, 
so  that  the  point  at  which  the  limit  cycle  velocity  goes  to 
zero  is  exactly  one-half  period  after  the  initial  point,  at 
-Xq.   In  this  case,  the  half-period  is  equal  to  the  constant 
T   and  T   equals  T  . 

Therefore,  the  linearization  in  which  we  are  interested 
is  completely  determined  by  the  first  half  of  the  limit 
cycle  trajectory,  because  of  the  odd  symmetry.   To  show 
this,  assume  that  the  same  deviation  from  the  point  X-,  is 
taken,  except  with  a  change  of  sign,  and  trace  the 
trajectory  from  this  point  x  -5x^.   Noting  that  this  initial 
point  equals  (-x  -fix  ) ,  and  using  the  odd  symmetry  of  the 
phase  portrait,  it  is  clear  that  the  trajectory  returns  to 
region  IV  at  the  point  (-x  -6x  )  (since  the  trajectory  from 
Xq+<5x   goes  to  x  +(Sx  )  .   Thus  the  resulting  deviation  on  the 
second  half  of  the  cycle  is 

(5.23)  (Sx^  =  (-Xj^-^Xj^)  -  Xq 

=  -% 
since  x.  =  -Xr,'   Linearizing  this  second  leg,  using  the 
initial  deviation  of  -5x  , 

(5.24)  <Sx  =  Z  (-<5x  )  +  higher  order  terms 

=  -% 

=  -Y  5x   +  h.o.t. 
by  the  analysis  of  the  first  half-cycle.   Since  this  holds 
true  for  any  deviation  Sx    ,    we  must  have  Y  =  Z.   Therefore, 

the  transition  matrix  for  the  whole  cycle  is  the  square  of 

2 
the  transition  matrix  for  the  first  half -cycle,  Y  . 
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Now  remove  the  assumption  of  symmetry.   Assume  that  the 
limit  cycle  is  not  symmetric  (but  still  involves  only 
sliding,  i.e.  case  1,  no  sticking);  in  this  case  the 
linearized  flow  map  must  be  examined  on  the  second  half  of 
the  cycle  also.   The  same  derivation  as  above  (equations 
(5.4)  through  (5.8))  can  be  performed,  but  starting  at  the 
initial  point  x,+5x  ,  in  order  to  determine  5x,^  in  terms  of 
Sx.  .   The  trajectory  in  region  II  is  found  using  autonomous 
state  variables: 

(5.25)  x  (t)  =  x(t)  +  1 

a 

so  the  trajectory  is 

(5.26)  x(t)  =  e^^^~'^l^  [x^  +  5x^  +  1]  -  1 

For  a  sufficiently  small  deviation  <5x  ,  the  velocity  goes  to 
zero  when  the  trajectory  reaches 

(5.27)  x(T  +  6T^)  =  e^^'^3'^'^'^3^  [x^  +  Sx^  +  1]  -  i 
The  deviation  of  this  point  from  the  point  where  the 
original  limit  cycle  reaches  region  IV  is 

(5.28)  <Sx^  =  x(T+5T3)  -  x^ 

=  {e^^T3  -  1}  e^^B  [x^  +  1]  +  e^^^S-^^^B^^x^ 

Note  that  by  comparing  this  formula  to  (5.8),  and  stating 
(5,8)  as  a  function  with  parameters 

(5.29)  6x^   =  {e^^T  -  1}  e^^l  [x^  -  1]  +  e^^V^^^^x^ 


f(5xQ,  ST;    T^,  Xq-I) 


then  we  see 


(5.30)     <S>Lj,   =  f(5x^,  6T3;  T3,  x^+1) 

The  linearization  for  this  function  has  already  been 

derived,  so  the  linear  map  for  the  second  half  of  the  cycle 
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is  found  by  substituting  the  correct  parameters  from  (5.3  0) 
into  (5. 19) ,  so  that 

(5.31)  Sx^   =  ^^^3'  ^1+1)  "^x^  +  h.  o.  t. 

(5.32)  Y(T3,  x^+1)   =  d(52Lp)/d(6x^)  |  ^^^^ 

=  e^'^3  -  (A  e^'^3  [x^+1]  }  (^'A  e^'^3  [x^+l]}""'-  {j0'e^'^3) 
By  the  argument  given  previously  in  (5.22),  the  total 
linearized  flow  map  is  the  product 

(5.33)  Y^^™^^^  =  Y(T_,  X.+l)  Y(T,  ,  x^-1) 

compos       3    1        1   ~0  ~ 

Sx^   =   Y        (Sx-  +  h.  o.  t. 
— T     compos   —0 

Note  that  this  reduces  to  the  symmetric  case  when  T   =  T 
and  x   =  -  x^: 

(^•34)     ^compos  =  ^(^1'  -^0+i)  ^(T^'  ^0-i) 

=  y2(T^,  Xq-1) 

since  Y  is  an  even  function  of  its  second  variable. 

The  final  step  in  the  derivation  for  this  case  is  to 

convert  back  to  physical  variables,  while  at  the  same  time 

forcing  the  initial  deviation  to  the  desired  region.   By 

following  the  reasoning  given  at  the  start  of  the  analysis 

for  this  case  (equation  (5.2)),  the  desired  matrix  for  the 

linearized  flow  map  of  the  entire  cycle  is 

(^•^^)     Vcompos  =  ^r(T3'  ^1^^)  ^r^^l'  ^0"^) 

Y^(T,,  Xq-I)  =  I^  T^2  ^(^1'  ^0-i)  ^12"'  ^r' 

(5.36)     6y^^  =   Y^_^^^p^^  ^y^^  +  higher  order  terms 
where  5y^   is  the  (n-1) -vector  representing  the  deviation  in 
physical  variables  after  going  around  one  cycle  (with 
velocity  =  Sy      left  out  since  it  is  zero) ,  5y    is  the 
initial  deviation  in  the  same  variables  (also  has  (n-1) 
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elements,  no  deviation  in  initial  velocity  allowed) ,  and  Y  , 

^r-compos  ^^®  ^^®  matrices  representing  the  linearized  flow 
maps  derived  above,  but  in  physical  variables  and  reduced  to 
order  (n-1) .   Note  that  the  latter  matrix  is  (n-l)x(n-l),  so 
its  (n-1)  eigenvalues  are  those  desired  for  the  local 
stability  analysis  (as  discussed  above,  the  n^'^  eigenvalue 
is  always  0) . 
Case  2;  Limit  Cycle  with  Sticking 

To  analyze  this  case,  it  is  necessary  to  derive  the 
linearized  flow  map  for  the  portions  of  the  limit  cycle  in 
the  sticky  region,  i.e.,  region  III.   These  linear  maps  can 
then  be  combined  with  the  maps  previously  derived  for  the 
sliding  regions  to  form  the  composite  map  for  the  cycle. 

Assume  an  initial  point  x^  that  meets  the  equations  in 
Chapter  IV,  so  that  this  is  the  point  at  which  the  system 
comes  to  rest  in  region  III  (this  is  the  fourth  switch  point 
in  the  cycle).   In  this  case,  the  period  the  system  sticks 
in  region  III  is  equal  to  the  constant  T   in  Chapter  IV. 
This  initial  point  must  be  used,  instead  of  the  point  x   at 
the  point  of  impending  motion,  in  order  to  derive  the 
behavior  within  region  III. 

A  deviation  Sx^    is  assumed  from  the  initial  point  of 
the  limit  cycle.   This  deviation  is  of  course  constrained  as 
before,  so  that  the  resulting  initial  point  remains  inside 
region  III  (velocity  =  0) .   Starting  at  the  initial  point 
^3''"*^3/  the  trajectory  is  determined  by  the  piecewise  linear 
methods  as  in  chapter  IV. 
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The  trajectory  in  region  III  is 

(5.37)  x(t)  =  e^s^  x(0)  =  e^s^  [X3  +  6X3] 

For  a  sufficiently  small  deviation  5x    ,    the  trajectory 
reaches  the  condition  of  impending  motion  (acceleration 
torque  sufficient  to  break  loose) ,  but  not  necessarily  in 
exactly  T   seconds.   Denoting  the  change  in  this  period  is 
<ST,  impending  motion  is  reached  when 

(5.38)  x(T^  +  ST)    =  e^s^^4"^*^^^  [x^  +  6X3] 

The  quantity  of  interest  is  the  deviation  of  this  point  from 
the  point  where  the  original  limit  cycle  breaks  free 

(5.39)  x^  =  e^s'^4  X3 

(5.40)  Sx^    =   x(T^+5T)  -  x^ 

,  A  (5T    ^,   A  T^     ^    A  (T^  +  5T)  ^ 
=  (es    -I}es4x   +   es^4     "^-i 

The  differentiation  of  this  function  to  determine  the 

first-order  dependence  of  5x   on  6x      is  performed  as  before; 

many  of  the  steps  are  similar  and  thus  omitted.   The  desired 

derivative  is 

(5.41)  Y^  =  d(5x^)/d(<5x3)  I^^^Q 

where  Y   is  the  desired  linearized  map  in  the  stiction 

region.   Note  that  since  x^  is  not  a  function  of  Sx^    or  ST, 

~4  —3 

(5.42)  [d/d(6x3)]  {Sx^}    =  [d/d(5x3)]  {x(T^+6T)  -  X4 } 

=  [d/d{Sx^) ]  (x(T^+5T) ) 

=  [d/d(6x3)]  (e^s(^4^'^^^  [X3  +  5X3]} 

The  dependence  of  Sx      on  both  <ST  and  Sx     explicitly  is 

handled  using  the  chain  rule,  as  before: 

(5.43)  d(5x^)/d(6x3)  =  D(5x^)/D(<!;x3) 

+  [D(6x^)/D(5T)]  [d(6T)/d(<SX3)] 
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where  "D"  is  again  the  symbol  for  partial  derivative.   The 
first  term  is 

(5.44)  D{Sx^)/D{6x^)       =    e^s'^4 

when  evaluated  at  Sx      =  ST  =   0.   The  partial  derivative  in 
the  second  term  is 

A     <?T 

(5.45)  D(5x^)/D(5T)       =    [D/D(<ST)]     {e    s         m) 

=  Ag  e  s  4  [X3] 
where  m  is  a  temporary  vector  independent  of  ST,    as  in  the 
previous  derivation. 

The  evaluation  of  the  [d  (<ST)/d  (<5x  )  ]  factor  in  the 
second  term  is  handled  as  before  by  relating  the  two 
variables  with  an  implicit  function,  and  using  the  implicit 
function  theorem.   The  implicit  relation  in  this  case  is 
determined  by  the  breakaway  condition  (i.e.,  impending 
motion) ;  as  in  the  analysis  in  Chapter  IV,  this  can  be 
stated  as 

(5.46)  §_'    A  x(T^+5T)  =  -Lg 

or,  equivalently, 

^'  A  e^s^'^4'^"^^^  [x^  +  5X3]  +  Lg  =  0 

Thus  we  have  a  function  of  the  form  f{ST,    Sx    )    =  0.   The 

derivative  is  then 

(5.47)  d(ST)/d{Sx^)    =      -    [Df/D(5T)  l"-*-  [Df/DiSx^)] 
Evaluating  these  factors, 

A  ST 

(5.48)  Df/D(5T)  =  [D/D(5T)]  {!'  A  e  s    m} 

o  I  TV  A    A  (ST 
=  fl'  A  A   e  s    m 
s       — 
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The  second  factor  is 
(5.49)     Df/D{Sx^)    =  [0/0(6X3)]  ^^'  ^  e^s^'^4"^"^'^^  5  X3 } 


=  £'  A  e  s^  4    ' 


so  that 


(5.50)  d(ST)/d(Sx^)    =  -{^'A  A^  e^s'^4  X3}"''"  {^'A  e^s'^4 } 
where  the  factors  are  again  evaluated  at  Sx=6T=0 .      Finally, 
combining  (5.44,45,50)  into  (5.43),  we  get  the  total 
derivative  representing  the  linearized  periodic  map 

(5.51)  Yg(T^,  X3)  =  d(5x^)/d(5x3)  | ^^^^ 

=  e^s'^4  -  (Ag  e^s'^4  X3}  {^'AAg  e^s'^4  X3}~''"  (^'A  e^s'^4 } 
Note  that  the  factor  to  be  inverted  is  a  scalar.   This 

expression  can  be  simplified  by  noting  that 

A  T 

(5.52)  X4  =  e  s  4  x_ 

so  substituting  into  (5.51)  gives 

(5.53)  Y^(T^,X3)  =  e^s^4  -  {A^  x^){£'AAg  x^}"-^(^'A  e^s'^4 ) 
=  e^s'^4  -  {Ag  x^  ^'A  e'^s^4)  /  {^'A  A^  x^ ) 

In  order  to  obtain  the  linear  map  for  the  entire  cycle 
for  the  case  with  sticking  during  the  cycle,  we  can  use  the 
chain  rule  as  before: 

(5-5^)     ^compos  =  (d(5x3)/d(5x2))  (d ( 5x2)/d ( 6x^) } 

X  (d(6x^)/d(5xQ) }  {d(5x^)/d(5x3)) 
where  the  last  factor  was  just  derived,  the  first  and  third 
are  obtained  from  the  sliding  case  previously  analyzed,  and 
the  second  factor  can  be  found  from  the  fourth  by  symmetry, 
as  shown  below.   Note  that  Sx      =  <5x  ,  since  both  are  the 
deviations  at  the  point  of  impending  motion  in  region  III; 
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this  equivalence  justifies  the  juxtaposition  of  the  last  two 
factors  when  the  chain  rule  is  applied. 

Note  also  that  in  (5.54)  the  order  of  the  factors  is 
different  from  that  used  in  the  sliding  friction  case;  the 
initial  deviation  for  the  case  with  sticking  is  taken  as  Sx 
at  the  point  x  where  the  system  comes  to  rest,  rather  than 
at  the  point  x  where  motion  is  impending.   This  order  is 
required  since  any  (small)  deviations  Sx     are  allowed  that 
keep  the  initial  velocity  zero  (hence  Sx      can  be  constructed 
as  in  (5.2)).   On  the  other  hand,  certain  initial 
perturbations  would  not  be  allowed  at  the  point  x  ;  any 
perturbations  at  this  point  of  impending  motion  that  altered 
the  forcing  torque  would  either  cause  motion  (so  initial 
point  no  longer  in  region  III) ,  or  would  mean  motion  is  no 
longer  impending.   This  problem  does  not  arise  in  the 
sliding  friction  case,  but  in  the  sticking  case  there  are 
two  constraints  on  Xq :  zero  velocity  and  forcing  torque  at 
breakaway.   It  is  therefore  much  more  convenient  to  take  the 
initial  point  to  be  x^  for  the  stability  analysis;  all 
constraints  are  then  taken  care  of  in  the  derivation 
automatically. 

The  second  factor  in  (5.54)  can  be  found  by  symmetry, 
as  follows:  taking  a  deviation  from  the  point  x   at  which 
the  system  comes  to  rest  in  region  III  of  Sx    ,    we  can  derive 
by  the  same  steps  as  above  that 
(5.55)     Sx^    =   x{T^+ST)    -    x^ 

=    (e^s^T  -  1}  e^sT2  x^  +   e^  (^2'^'^^)  5x^ 
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The  similarity  of  this  equation  to  (5.40)  can  be  used  to 
show  that  the  desired  linearization  is 
(5.56)     Y^(T2,  x^)  =   d(Sx^)/d{Sx^)     \ ^^^^ 

=   e^s'^2  -  {Ag  e^s^2  x^ }  {^'AA^  e^s'^2  x^)"^  {^'A  e^s'^2 } 
Of  course,  as  in  the  sliding  case,  this  matrix  must  be 
modified  to  transform  to  physical  variables,  and  reduced  to 
order  (n-1) . 

The  linear  transformation  for  the  entire,  four-part 
cycle  is  then 

^^■^'^  Vcompos  =  Yr(T3'^2-^i)  ^rs^^s'^) 

X  Y^(T^,Xo-l)  ^rs^T^^x^) 


Y 


-1 


r(T,,  Xq-1)  =  I^  T^2  ^(^1'  ^0-i)  T^2    ^r 
Y^3(T2,  X,)  =  I^  T^2  Y^CT^,  X,)  T^^""  ^r' 
(5.58)     Sy^^  =   Y^_^^^p^^  Sy^^   +  higher  order  terms 
where,  as  before,  the  5y^  vectors  are  the  (n-1)  elements  of 

the  initial  and  final  deviation,  respectively,  and  Y 

-^        r-compos 

is  the  desired  linear  transformation  that  determines  local 
stability  of  the  limit  cycle. 

The  results  of  the  stability  analysis  of  this  chapter 
can  be  summarized  as  follows: 

Theorem  5-1;   Given  the  existence  of  a  friction  limit  cycle, 
as  predicted  by  the  solution  of  the  nonlinear  algebraic 
equations  derived  in  Chapter  IV  (Theorem  4-1) ,  the  local 
asymptotic  stability  of  the  periodic  orbit  is  determined  by 
the  eigenvalues  of  the  matrix: 
(Case  1:  Sliding  only) 

(^•^^^     Vcompos  =  ^r(T3'  ^1+^)  ^r^^l'  ^o'^^ 
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(Case  2:  Sliding  plus  sticking) 

X  Y^(T^,  Xq-1)  y^^(T^,  X3) 


where 


Y^(T,,  Xq-1)  =  I^  T^2  ^(^1'  ^0-i)  ^12""  ^:  ' 


r 
r 


Y^3(T2,  x^)  =  I^  T^2  ^s^T2'  ^1)  ^12    ^. 
the  x^'s  are  the  switching  points,  the  T.'s  are  the 
corresponding  switching  periods,  I   is  the  (n-1)  x  n  matrix 
defined  in  (5.3),  T^^  i^  the  transformation  from  physical  to 
canonical  variables  defined  in  Chapter  III,  and  the  matrix 

functions  Y  and  Y   are 

s 

(5.19)     Y(T^,  Xq-1)  =  d(5x^)/d(5Xo)  I 5x-0 

=   e^*^!   -    {A   e-^^l    [Xq-I]  }     {jS'A   e-^"^!    [Xq-1  ]  )  "^{I '  e^^l } 

(5.56)  Y^(T2,    x^)     =    d(5x2)/d(5x^)      | ^^^^ 

AT  AT  a    T  —1  a    t> 

=    e    s    2    -    {A^    e"s^2    x^ }     {^'AA^    e    s^2    x^ }    ^    {^'A   e^s^2} 

Before  this  theorem  is  proved,  some  terms  must  be 
defined  more  rigorously. 

Definition:   The  flow  of  the  differential  equation  system  is 
the  map  (t , x^) : $ (t , x^) ,  where  i(t,x^)    =   x(t)  when  the 
initial  condition  x(0)  =  x  . 

Definition;   The  Poincare  map  (P-map)  is  the  map  6x  :g(<5x  ), 
where  g(<5xQ)  =  x(ti)  -  x^,  5x^   =   x(0)  -  x^ ,  and  tj  is  the 
first  time  at  which  the  trajectory  again  crosses  the  local 
section  around  the  periodic  orbit  at  x  . 

Proof:  The  proof  follows  from  the  theory  in  Chapter  13  of 
Hirsch  and  Smale  (1974)  and  the  calculations  above,  except 
for  one  difficulty  to  be  overcome.   The  reference  theory  is 
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stated  for  the  case  where  f(.)  is  continuously 
dif ferentiable  (so  that  the  flow  *  is  continuous) .   However, 
as  will  be  shown  below,  all  that  is  required  to  obtain  the 
desired  result  is  that  the  flow  be  a  continuous  function  of 
the  deviation  about  the  initial  point,  x  +  Sx    ,    so  that 
certain  limits  can  be  taken,  and,  in  addition,  that  the 
P-map  is  dif ferentiable,  so  that  the  derivative  may  be  used 
to  check  asymptotic  stability. 

As  a  preliminary  step,  then,  these  properties  must  be 
shown.   The  flow  is  the  integral  of  a  real-valued  function 
of  time  (see  equation  5.69  below),  f(x(,)),  so  it  is  a 
continuous  function  of  time  (Rudin,  (1976),  Theorem  6.20). 
It  is  also  clear  that  the  flow  is  continuous  in  x  within 
each  piecewise  linear  region,  since  the  function  f(.)  is  at 
least  C-'-  there  (and  Hirsch  and  Smale,  Chapter  15,  shows  this 
is  sufficient  for  continuity  of  the  flow) .   Now,  taking  the 
initial  condition  as  a  sufficiently  small  deviation  from  the 
initial  point  x  ,  the  flow  at  any  finite  time  is  the 
composition  of  a  finite  number  of  continuous  functions, 
formed  from  the  pieces  of  the  trajectory  in  each  region  it 
traverses.   This  composition  can  be  formed  since  the 
trajectory  is  continuous  across  the  discontinuities  in  f ( . ) ; 
that  is,  xCt"*")  =  x(t~)  ,  where  t  is  a  switching  time. 
Therefore  the  flow  is  a  continuous  function  of  this  initial 
condition  (Rudin,  (1976),  Theorem  4.7,  the  composition  of 
continuous  functions  is  continuous) ,  for  sufficiently  small 
deviations. 
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The  same  approach  is  used  to  demonstrate  the 

differentiability  of  the  P-map.   This  map  is  simply  the 

composition  of  the  flows  in  each  of  the  regions,  so  is 

continuous  by  the  previous  argument.   But,  in  fact,  the 

flows  in  each  region  are  continuously  dif ferentiable,  so 

their  composition  is  also  dif ferentiable. 

Using  these  properties,  the  rest  of  the  proof  proceeds 

as  in  Hirsch  and  Smale  (1974),   In  chapter  13,  section  3,  it 

is  shown  that  if  the  P-map  has  x   as  a  sink,  then  the  orbit 

is  asymptotically  stable  (using  the  property  of  continuity 

of  the  flow  to  show  that  the  flow  converges  uniformly  to  the 

periodic  orbit  if  the  deviations  after  each  cycle  also 

converge) .   But  the  P-map  is  a  discrete-time  dynamical 

system,  so  it  is  a  sink  if  the  derivative  of  this  map  has 

all  eigenvalues  less  than  one  in  magnitude  (which  requires 

the  property  of  differentiability  of  the  P-map) .   The 

calculations  of  this  derivative  have  been  presented  above, 

and  show  that  the  matrix  Y         is  in  fact  the  derivative 

r-compos 

of  the  P-map  at  x  .   Therefore  the  (n-1)  eigenvalues  of  this 
matrix  determine  if  the  orbit  is  asymptotically  stable,  and 
the  proof  is  complete. 

To  summarize  the  above  stability  results,  it  has  been 
shown  that  the  perturbations  from  the  periodic  orbit  evolve 
according  to  the  equation 

(5.58)     5:^^  =   Y^_^^^p^^  SYq^      +  higher  order  terms 
where  ^y,^   is  the  (n-1) -vector  representing  the  deviation  in 
physical  variables  after  going  around  one  cycle  (with 
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velocity  =  5y^  left  out  since  it  is  zero)  ,  and  5y_        is  the 

initial  deviation  in  the  same  variables  (also  has  (n-1) 

elements,  no  deviation  in  initial  velocity  allowed) .   The 

eigenvalues  of  the  matrix  Y        ,  therefore,  indicate  the 

r-compos 

stability  of  the  periodic  orbit  (i.e.,  if  they  are  all  less 
than  one  in  magnitude,  the  orbit  is  asymptotically  stable, 
and  in  fact  a  periodic  attractor) . 

In  addition,  it  has  been  shown  that  the  results  on 
orbital  stability  in  Hirsch  and  Smale  (1974),  chapter  13, 
can  be  extended  to  the  case  where  f(.)  is  only  piecewise 
linear,  rather  than  C^. 

Stability  Calculations  for  Example  Problems 

The  limit  cycle  orbital  stability  equations  derived 
above  can  be  applied  to  the  example  system  discussed  in 
Chapters  I  and  IV  (example  1-2,  case  2,  a  3D  system  with 
both  a  sliding  and  stiction  limit  cycle) .   The  example 
system  is  defined,  as  in  the  example,  by 
(5.59)  dy/dt   =   A  y  +  b  L^ 


where 


A   = 


0 
0 
-K, 


1 

0 

-K, 


0 
1 
-B 


1      2 
with  b'  =  [0  0  -1],  L^  =  friction  torque,  and  where  the 

state  variables  are  y   =  position,  y  =  velocity,  and 

y^  =  compensator  integrator.   The  transformation  matrices  T 

and  T   are  identity  matrices,  and  the  row  vector  ^'  equals 

[0  0  1].   The  particular  case  has  B  =  0.9,  K.  =  K  =  1,  with 
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closed-loop  poles  at  .026  ±  j  1.024  and  -.9524.   The 
stiction  limit  cycle  solution  is  at  approximately 
X   =  [0.98  0.22  0],  with  T   at  approximately  n .    There  is 
also  a  solution  of  the  equations  for  the  symmetric,  sliding 


case,  at  T   =  3.1455  seconds,  and  x   =  [-0.283  12.4  0]. 
Example  V-1,  Part  1:   Stability  of  Sliding  Limit  Cycle 

For  this  case,  equations  (5.35,  5.19)  from  Theorem  5-1 
apply.   Using  the  parameters  defined  above,  we  find  that 


(5.60) 
so  that 


AT 
Ae   1  [X   -  1]  =  [-12.41  .001  13.13]' 


AT 

§_*    Ae   1  [X   -  1]  =  13.13 


The  matrix  exponential  is  (approximately,  a  numerical 


algorithm  was  used) 
AT. 


(5.61) 


1   = 


-.5314 
-.5227 
.5830 


-.1126 

-1.054 
.0603 


.5227 

-.5830 

-.5294 


Therefore,  the  matrix  Y  required  to  define  the  linearized 
flow  map  is 

(5.62)     Y   =  .0197    -.0556     .0223 

-.5227   -1.054     -.583 
.0002    .0000      -.0002 
Although  there  are  unavoidable  numerical  errors  in  this 
calculation,  the  true  Y  matrix  apparently  does  have  a  zero 
eigenvalue,  as  expected.   The  upper  left  block  sub-matrix  is 
equal  to  Y  ,  and  the  composite  stability  matrix  is 


(5.63) 


02944  .05749 
5407   1.14 
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This  matrix  has  eigenvalues  of  .002125  (with  corresponding 
eigenvector  at  [-2.104  1.]'  ),  and  1.168  (with  eigenvector 
[.05  1.]'  ).   Thus  it  has  one  stable  and  one  unstable  mode 
(note  this  is  a  discrete-time  system) ,  and  is  a  saddle 
point.   A  simulation  of  this  system  was  used  to  verify  this, 
where  the  x   and  x_  states  were  sampled  whenever  the 
velocity  x  went  through  zero  (near  x_) .   The  plot  in  Figure 
V-1  shows  the  point  x  ,  the  eigenvectors  of  the  stable  and 
unstable  modes,  and  typical  trajectories  near  x  .   The 
stable  mode  is  very  fast  so  that  simulated  trajectories 
quickly  jump  onto  the  unstable  eigenvector,  then  move  along 
it.   The  simulation  experimentally  verifies  the  stability 
predictions. 
Example  V-1,  Part  2;   Stability  of  Sticking  Limit  Cycle 

Similar  calculations  can  be  performed  for  the  other 
limit  cycle  that  exists  for  this  system,  except  that 
equation  (5,57)  from  Theorem  5-1  governs  this  case.   Using 
the  initial  condition  defined  above,  with  T   =  tt  and  T   =  9 
(as  approximate  values),  and  x  =  x(T  )  =  [1.  -0.22   0.]', 
and  recalling  that  the  symmetry  defines  the  other  switching 
points  and  periods,  we  find  after  substitution  in  (5.56) 

0      -1    -9.9 

0      19 

0      0     1 
with  Y   being  the  upper  left  block  sub-matrix.   The 
stability  matrix  for  the  sliding  portions  of  the  limit  cycle 
is  also  required,  and  the  upper  left  block  of  Y  is 


(5.64)     Y 
^     '      s 
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(5.65)     Y^    =        .02034   -.05551 

-.5267   -1.055 

The  composite  stability  matrix  representing  the  linearized 

flow  map  for  the  whole  cycle  is 

2 


(5.66) 


(  y   Y    ) 
^   r  rs  ' 


0      .0400 

0      .2786 

which  has  eigenvalues  0.  (eigenvector  [1  0]')  and  .2786 
(eigenvector  [.144  1.]').   This  limit  cycle  is  stable. 
Simulation  results,  obtained  as  above,  are  shown  in  Figure 
V-2 .   Note  that  the  zero  eigenvalue  causes  the  trajectories 
to  jump  to  the  eigenvector  for  the  non-zero  eigenvalue, 
where  they  then  converge  to  x  .   These  results  match  the 
predictions  of  the  stability  equations  for  this  stiction 
case,  also. 

The  global  phase  protrait  can  now  be  qualitatively 
inferred  from  the  known  limit  cycle  orbits  and  their 
stability,  if  it  is  assumed  no  other  periodic  solutions 
exist.   Simulation  results  are  displayed  in  Figure  V-3 
through  V-5,  confirming  this  picture.   Trajectories  outside 
the  unstable,  sliding  mode,  limit  cycle  (that  is,  outside 
the  stable  manifold  of  the  unstable  I.e.)  will  spiral  to 
infinity.   This  stable  manifold  therefore  represents  a 
stability  boundary,  since  trajectories  inside  it  will  never 
leave  it  (it  is  an  invariant  surface) .   There  is  a  line 
segment,  including  zero,  that  are  equilibrium  points,  since 
for  x   =  X  =  0  and  for  x   less  than  Ls  in  magnitude,  there 
is  insufficient  torque  to  breakaway  (they  are  unstable  since 
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any  disturbance  in  x^  or  x^  would  cause  them  to  spiral 
away) .   Trajectories  starting  at  any  points  inside  the 
stability  boundary,  except  for  on  this  line  segment,  will 
approach  the  stiction  limit  cycle  asymptotically. 
Symmetry  of  Limit  Cycles  for  Odd  Systems 
This  section  discusses  the  question  of  the  symmetry  of 
any  existing  limit  cycles  of  a  dynamical  system  with  an  odd 
time  deriyative.   That  is,  if  the  system  is  described  by 

(5.67)  dx/dt  =  f(x) 

and  the  function  f(.)  is  odd,  are  all  the  limit  cycles  of 
the  system  necessarily  half-waye  symmetric, 

(5.68)  x(t  +  T/2)  =  -x(t) 

for  all  t,  with  T  the  limit  cycle  period?   The  answer  to 
this  question  is  no,  except  possibly  for  certain  specific 
cases,  as  will  be  seen  below,  but  it  is  interesting  to 
examine  what  conditions  are  required  to  proye  symmetry,  and 
other  related  issues. 

These  issues  are  significant,  since  most  nonlinearities 
of  interest  in  seryo  design  are  odd,  piecewise  linear 
functions  (backlash,  friction,  quantization,  limiters, 
saturation,  deadzone,  .  .  .),   Also,  this  greatly  simplifies 
the  task  of  finding  friction  limit  cycles  by  the  piecewise 
linear  method  (the  topic  of  this  thesis) ,  since  it  is  only 
necessary  to  calculate  half  the  trajectory  of  any  limit 
cycle.   As  the  situation  stands  now,  if  a  trajectory  is 
found  that  maps  some  initial  point  x  into  -x,  it  is  known 
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that  a  limit  cycle  exists  with  twice  that  period,  but  the 
converse  is  not  necessarily  true. 

Several  results  about  this  problem  are  presented, 
although  the  general  problem  is  still  open.   Several 
counter-examples  are  given  for  various  specific  conjectures 
about  the  relationship  between  oddness  and  symmetry,  but 
symmetry  is  proved  for  systems  of  second  order. 
N-dimensional  results  are  also  given,  as  are  some  specific 
properties  of  friction  limit  cycles. 
Results  for  General  Odd  Systems 

The  condition  (5.68)  is  equivalent  to  the  statement 
that  for  every  x  belonging  to  the  trajectory,  -x  also 
belongs  to  the  same  trajectory  (that  is,  the  cycle  is 
symmetric  about  the  origin) ,  as  can  be  seen  by  looking  at 
the  solution  to  (5.67): 

(5.69)  x(t)  =  x(0)  + 
and  using  the  oddness  of  f(.).   Thus,  the  time  symmetry  of 

(5.68)  is  the  same  property  as  the  geometrical  symmetry  of 
the  periodic  orbit  about  the  origin. 

If  the  nonlinearity  is  odd,  some  information  about  the 
structure  of  the  phase  portrait  is  known: 
Theorem  5-2:   If  f  is  odd  in  (5.67),  and  x(t)  is  a 
trajectory  of  (5.67),  then  so  is  the  curve  -x(t) . 
Proof:   Since  x(t)  is  a  trajectory,  we  have  by  (5.69) 

(5.70)  x(t)  =  x(0)  +   Jq  f(x(s))  ds 
for  some  initial  condition  x(0) .   Look  at 

(5.71)  -x(t)  =  -x(0)  -  Jq  f(x(s))  ds 


nt 
Q  f(x(s))  ds 
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=  -x(0)  +  Jq  f(-x(s))  ds 
where  the  oddness  of  f(.)  is  used  in  the  last  step. 
Equation  (5.71)  shows  the  trajectory  -x(t)  is  a  solution  for 
(5.67),  starting  at  -x(0);  but  the  solution  defines  the 
trajectory  uniquely  and  the  proof  is  complete. 

Therefore,  it  can  be  concluded  that  the  phase  portrait 
is  symmetric  about  the  origin,  so  that  any  half-space  can  be 
found  from  the  other  half-space  by  reflection.   Although 
this  seems  very  close  to  the  desired  property,  in  fact  there 
are  significant  problems  in  proving  half-wave  symmetry  (or 
symmetry  of  limit  cycle  about  origin)  from  phase  portrait 
symmetry  and  the  oddness  of  f(.).   Several  authors  have 
attempted  to  use  this  fact  to  prove  that  the  limit  cycles  of 
a  system  with  an  odd  nonlinearity  are  symmetric  with  respect 
to  the  origin  (for  example,  Kovatch  [1]).   A  defect  of  these 
proofs  is  that  it  is  not  clearly  stated  that  the  result 
depends  on  the  specific  structure  of  the  nonlinearities 
assumed  by  those  researchers.   In  fact,  oddness  is  not 
enough  to  prove  symmetrical  limit  cycles,  as  shown  by  the 
following  counterexample: 
Example  V-2 :   A  system  has  the  form 
(5.72)     dx^/dt  =  X2 

dx  /dt  =  sin(x  ) 
Then  the  Jacobian  of  the  system  at  (27r,  0)  is 
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(5.73) 


-cos(x  ) 


1 
0 


0 

-1 


1 
0 


whose  eigenvalues  are  j , - j ,  which  implies  that  this  is  a 
center.   In  fact  there  are  an  infinite  number  of  periodic 
solutions  near  this  point,  which  are  not  symmetric  with 
respect  to  the  origin.   Yet  the  system  function  f(.)  is  odd! 

This  2D  system  is  odd  and  has  limit  cycles  which  are 
not  half-wave  symmetric.   The  phase  portrait  is  symmetric 
about  the  origin  as  required  by  the  theorem  described  above, 
but  there  is  not  a  unique  limit  cycle.   This  leads  to: 
Theorem  5-3;   If  the  limit  cycle  for  (5.67)  is  unique  with 
f(.)  odd,  then  the  limit  cycle  is  half -wave  symmetric. 
Proof:   We  know  by  previous  theorem  that  phase  portrait  is 
symmetric.  Therefore,  if  x(t)  is  a  limit  cycle,  -x(t)  is 
also  a  closed  periodic  trajectory.   If  x(t)  is  not  symmetric 
about  the  origin  (equivalent  to  half-wave  symmetric)  then 
reflecting  about  the  origin  generates  a  separate  limit 
cycle,  contradicting  the  hypothesis.   This  completes  the 
proof . 

In  addition,  by  a  simple  extension  of  the  above  proof, 
one  can  show: 

Theorem  5-4:   If  the  limit  cycle  for  (5.67)  is  unique  in  any 
region  that  is  symmetric  about  the  origin,  with  f(.)  odd, 
then  the  limit  cycle  is  half-wave  symmetric. 
Proof;   By  a  similar  argument,  reflecting  the  limit  cycle 
about  the  origin  generates  another  limit  cycle  which  must  be 
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contained  in  the  region  (since  it  is  symmetric  also) ,  which 
leads  to  a  contradiction. 

Thus  a  limit  cycle  known  to  be  unique  in  a  sphere  or 
spherical  shell,  for  example,  must  be  symmetric. 

Further  examination  of  Example  V-2  shows  that  in  fact 
the  limit  cycle  is  symmetric  about  its  associated 
equilibrium  point,  which  leads  to  the  conjecture  that  if 
f(.)  is  odd  in  (5.67),  then  any  limit  cycles  are  symmetric 
with  respect  to  an  equilibrium  point.   Equivalently ,  making 
a  change  of  variable  so  that  the  equilibrium  point  is  at  the 
origin  (as  long  as  f  is  still  odd) ,  then  the  limit  cycle  is 
half-wave  symmetric,   A  slightly  different  result  that  might 
be  easier  to  show  would  be 

Conjecture:   If  (5.67)  has  an  equilibrium  point  at  the 
origin,  f(.)  is  odd,  and  a  limit  cycle  exists,  then  the 
limit  cycle  is  symmetric. 

In  spite  of  Example  V-2,  the  conjecture  is  true  for  the 
special  case  of  2D,  with  the  extra  condition  of  a  unique 
equilibrium  point: 

Theorem  5-5:   If  (5.67)  is  a  2D  system  and  has  a  unique 
equilibrium  point  at  the  origin,  then  the  conjecture  holds. 
Proof:   A  limit  cycle  in  a  2D  dynamical  system  must  contain 
an  equilibrium  point  (using  index  theory,  index  of  closed 
curve  must  be  one;  see  Vidyasagar  (1978) ,  Theorem  69  in 
section  2.4).   By  the  hypothesis  in  the  conjecture,  the 
equilibrium  point  inside  the  limit  cycle  is  located  at  the 
origin.   Since  (5.67)  is  odd,  the  reflection  of  the  limit 
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cycle  exists  and  is  a  limit  cycle  containing  the  origin  (the 

origin  would  map  onto  itself  and  thus  be  in  the  interior  of 

both  limit  cycles) .   Then  taking  the  point  x   that  achieves 

s 

the 

(5.74)  sup  I  X  I  =  x^^ 
(X  e  LC) 

(that  is,  the  point  on  the  limit  cycle  that  is  at  the 
maximum  radius  from  the  origin)  and  the  point  x.  that 
achieves  the 

(5.75)  inf  I  X  I  =  X 
(X  e  LC) 

(which  both  exist  and  are  achieved  by  continuity) ,  it  is 

seen  that  the  reflected  limit  cycle  has  the  same  points,  but 

at  -x^  and  -x . .   Then  -x.  must  be  in  the  interior  or 
— s      —1         —1 

touching  the  trajectory  x(t) ,  and  -x   must  be  in  the 
exterior  or  touching,  x(t)  divides  the  plane  into  distinct 
regions,  and  the  trajectory  -x(t)  is  continuous.   Therefore, 
the  trajectories  must  touch  or  cross,  or  be  identical.   The 
latter  result  is  the  only  one  not  leading  to  a 
contradiction . 

This  result  could  probably  be  extended  to  cover  the 
case  of  Example  V-2  also,  if  the  condition  was  weakened  to 
require  only  unique  equilibrium  within  a  symmetric  region, 
for  example. 

Unfortunately,  a  counter-example  exists  for  the  general 
case  of  n-dimensions,  so  the  conjecture  does  not  hold. 
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Example  V-3:  (Found  by  V.  M.  Popov,  personal  communication) 
Consider  the  following  system: 

(5.76)  dx/dt  =   -(x^  +  y^  -  1)  X   -   y 
dy/dt   =   X   -   (x^  +  y^  -  1)  y 

dz/dt   =   -(z^  -  1)  z  -  (x^  +  y^  -  z^) (X  +  y  +  z) 
This  system  has  a  single  equilibrium  point,  the  origin, 
which  is  unstable  (eigenvalues  of  the  Jacobian  at  the  origin 
are  +1,  1-jl,  and  1+jl).   It  also  has  two  periodic 
solutions,  at 

(5.77)  x(t)  =  cos(t),   y(t)  =  sin(t),   z (t)  =  +1  or  -1 
Examination  of  the  differential  equations  show  that  the 
system  is  odd,  since  all  terms  on  the  right-hand  side  are 
first  or  third  order  in  x,  y,  and  z.   Therefore,  an  odd 
system  function  and  a  unique  equilibrium  point  is 
insufficient  to  guarantee  limit  cycle  symmetry  (with  respect 
to  the  origin)  for  n>2.   Note,  however,  that  the  phase 
portrait  is  symmetric,  since  one  limit  cycle  is  the 
reflection  of  the  other. 

Results  for  Lure-Tvpe  Systems 

Some  results  have  been  completed  for  the  case  of  a 
linear  system  in  combination  with  a  single  SISO  nonlinearity 
(called  a  Lure-type  system).   First  of  all,  if  the  linear 
portion  of  the  system  (that  is,  the  transfer  function 
between  the  nonlinearity  output  and  input,  forming  a 
feedback  loop)  is  described  by  a  transfer  function,  or,  a 
state  variable  representation  that  is  controllable,  it  can 
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be  realized  in  a  controller  canonical  form,  as  shown  in 
Chapter  III  for  friction  limit  cycles: 

(5.78)  dx^/dt  =  X2 
dx^/dt  =  X 

dx   ,/dt    =  X 
n-1^         n 

dx  /dt  =  a   Xt  +  . . .  +  a   .  x   +  nonlinear  term 
n       01  n-l   n 

The  development  of  the  transformation  to  control  form  in 

Chapter  III  thus  allows  the  following  result  to  be  obtained: 

Theorem  5-6:   All  equilibrium  points  and  the  centroids  of 

all  limit  cycles  are  on  the  line  x„  =  . . .  =  x   =0  (i.e., 

2  n      ^     ' 

the  X- -axis) . 

Proof:   Equilibrium  points  are  defined  by  the  points  where 
f^(x)  =  f2(x)  =  ...  =  f^i^-^  "  °'  ^^^  first  (n-l)  of  these 
requirements  prove  the  equilibrium  point  result.   If  a  limit 
cycle  exists,  we  have  (by  (5.69)  and  the  limit  cycle 
condition  that  x(T)  =  x(0) 

(5.79)  Xq  =  Xq   +   Jq  f(x(s))  ds 


Therefore 

T 

(5.80)  Jq  f(x(s))  ds  =  0 

But  the  i    component  of  the  centroid  is  defined  by 

T 

(5.81)  x^.  -  (1/T)   Jq  x.(s)  ds 

Substituting  f   through  f  _   from  (5.78)  into  (5.80)  shows 
(5.81)  is  zero  for  x   through  x  ,  and  the  proof  is  done. 

However,  controllability  of  the  linear  portion  of  the 
system  is  insufficient  to  guarantee  symmetry.   One  need  look 
no  farther  than  Example  V-2  to  find  a  counter-example:   this 
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system  is  in  the  control  canonical  form,  but  has  equilibrium 
points  and  centroids  of  periodic  cycles  distributed  along 
the  entire  x   -  axis!   Another  counter-example  can  be 
constructed  as  follows: 

Example  V-4 :   Given  any  system  with  a  symmetric  limit  cycle, 
a  system  can  be  generated  that  has  an  infinite  number  of 
limit  cycles,  merely  by  augmenting  the  state  model  with  a 
single  new  state  whose  derivative  is  the  state  x  . 
Renumbering  the  original  states  as  2  through  (n+1) ,  and 
letting  the  new  state  be  x  ,  the  system  (5.78)  becomes 
(5.82)     dx^/dt  =  x^ 
dx^/dt  =  x^ 


dx  /dt  =  X  , ^ 
n'^       n+1 


X  ,,  +  nonlinear  term 


The  system  is  still  in  the  controller  form,  and  therefore  is 
controllable  (in  the  linear  portion) .   The  system  still  has 
the  same  limit  cycle  since  the  new  state  is  simply  the  time 
integral  of  the  state  x   (old  state  x  ),  and  the 
zero-centroid  property  for  the  original  symmetric  limit 
cycle  guarantees  this  time  integral  is  periodic.   However, 
the  same  conditions  hold  for  any  initial  condition  of  the 
state  X  ;  the  limit  cycle  returns  to  this  initial  condition 
after  one  cycle.   Thus  an  uncountably  infinite  number  of 
identically-shaped  limit  cycles  exist,  forming  a 
(hyper-) cylinder  parallel  to  the  x  -axis!   This  example 
system  has  a  so-called  free  integrator,  and  hence  has  a  zero 
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eigenvalue  in  the  closed-loop  linear  system  (or  in  fact  in 
the  closed-loop  nonlinear  system) . 

Clearly,  some  kind  of  observability  condition  is  also 
required;  in  example  V-3  the  state  x   could  not  be  observed 
from  the  system  output  (the  input  to  the  nonlinear ity  for 
the  Lure  system) .   Therefore,  controllability  and 
observability  of  the  linear  portion  is  necessary  to 
guarantee  symmetry.   Since  the  system  in  Example  V-2  meets 
these  requirements,  it  can  be  seen  that  a  condition  on  the 
location  and/or  uniqueness  of  equilibrium  points  is  also 
necessary.   A  conjecture  about  the  symmetry  of  limit  cycles 
for  odd  systems  containing  these  requirements  can  now  be 
stated: 

Conjecture;  Given  a  Lure  type  nonlinear  system  with  an  odd 
nonlinearity,  any  existing  limit  cycles  would  be  symmetric, 
if  the  system  meets  the  following  conditions: 

(1)  A  unique  equilibrium  point  (or  possibly  line 
segment)  exists  at  the  origin,  and 

(2)  The  linear  portion  of  the  system  is  represented  by 
a  transfer  function,  or  a  controllable  and  observable  state 
variable  model. 

This  remains  an  important  result  if  proved,  despite  the 
restrictions,  since  nonlinearities  of  interest  in  the  design 
of  servomechanisms  (such  as  those  listed  at  the  beginning  of 
this  section)  generally  meet  these  requirements.   The 
conjecture  is  still  open  at  this  time. 
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Specific  Results  for  Friction  Limit  Cycles 

There  is  still  hope,  then,  that  it  can  be  proved  that 
friction  limit  cycles  are  always  symmetric,  despite  the 
failure  of  the  conjecture  for  the  general  case  of  odd 
systems.   This  nonlinearity  meets  the  conditions  of  the 
conjecture,  since  it  is  odd  and  the  equilibrium  points  are  a 
line  segment  containing  the  origin  along  the  x   axis. 

This  subsection  examines  the  specific  case  of  the 
friction  nonlinearity  for  symmetry  results.   The  equations 
previously  derived  for  friction  limit  cycles  can  be  analyzed 
for  information  about  friction  limit  cycle  symmetry;  these 
results  would  not  apply  for  general  odd  systems,  although 
they  may  be  extendable  beyond  the  friction  case. 
Theorem  5-7;   A  (purely  sliding)  friction  limit  cycle  is 
asymmetric  if  and  only  if  the  switching  periods  (T   and  T   n 
previous  chapters)  are  unequal.   Sufficiency  holds  for  the 

case  where  the  system  matrix  A  has  no  purely  imaginary  (or 

AT 
zero)  eigenvalues,  and  therefore  {I  -  e   }  is  invertible. 

Proof:   Suppose  a  limit  cycle  exists  that  is  asymmetric,  and 

meets  the  necessary  conditions  of  Chapter  IV  for  the  purely 

sliding  case.   By  the  property  of  odd  systems  that  the  phase 

portrait  is  symmetric,  there  must  exist  a  distinct  limit 

cycle  that  is  the  reflection  of  the  first.   By  the  equations 

in  Chapter  IV,  these  two  solutions  have  initial  conditions: 

(5.83)     Xq   =   g(T3,  T^) 

=  -{I  -  e^*^)"^  (e^*^  -  2e^'^3  +  I)  1 


Ill 


Yq  =  g(Ti,  T3) 

=  -{I  -  e^"^}"^  {e^"^  -  26^"^!  +  1}  1 
If  T  =  T  ,  the  initial  conditions  are  the  same  and  the 
limit  cycles  are  not  distinct,  a  contradiction.   This  proves 
necessity. 

For  the  other  direction,  assume  a  limit  cycle  with  T 
not  equal  T^ 7  it  would  have  the  initial  condition  x   defined 
by  the  equation  above.   Assume,  by  way  of  contradiction, 
that  the  limit  cycle  is  symmetric.   Then  the  reflected  limit 
cycle,  which  should  be  the  same  trajectory  by  symmetry, 
would  have  initial  condition  y   above.   Since  these  points 
should  be  identical,  we  have 
(5.84)     0   =   ^0  ~  ^0 

=   {I  -  e^'^}"-^  {2e^'^3  -  2e'^^l}  1 

0   =   {e^'^3  -  e^"^!}  1 

,  A5T    T,   AT,  , 
=   {e     -  I)  e   11 

The  matrix  exponential  cannot  have  a  zero  eigenvalue,  so 

AT     .  A  5  T 

e   1  1  IS  a  nonzero  vector.   Therefore  e    must  have  an 

eigenvalue  of  1.   But  6T  =  T   -  T   is  nonzero,  so  the 

corresponding  eigenvalue  of  A  must  have  a  zero  real  part, 

which  contradicts  the  hypothesis.   This  proves  sufficiency, 

completing  the  proof. 

Note  that  this  proof  leaves  open  the  possibility  that  a 

system  with  an  imaginary  pair  of  eigenvalues,  with  periods 

matching  the  limit  cycle  period,  may  have  asymmetric  limit 

cycles.   This  would  be  similar  to  the  case  of  example  V-3 ; 

an  example  might  be  formed  by  augmenting  a  system  with  a 
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limit  cycle  appropriately.   Such  a  system  would  have 
multiple  (possibly  an  infinite  number  of)  limit  cycles 
related  to  the  original  one,  but  distinct  from  it  by 
oscillatory  terms  in  these  extra  states  (examples  exist  of 
systems  with  limit  cycles  that  are  entwined  about  each 
other,  such  as  the  Duffing  equation) .   The  example  would  be 
highly  contrived,  however,  since  it  would  have  to  meet  the 
controllability/ observability  conditions  with  the  new 
augmented  system,  yet  not  affect  the  existence  of  the  limit 
cycle.   No  such  example  has  yet  been  constructed.   This  case 
seems  somewhat  irrelevant  to  the  problems  typically  found  in 
servo  design;  therefore,  this  case  can  be  eliminated  by 
assuming  no  eigenvalues  with  zero  real  part.   Note  that  this 
also  eliminates  the  free  integrator  case,  and  is  consistent 
with  assumptions  made  in  Chapter  III. 

Combining  this  theorem  with  a  previous  theorem  on  the 
centroids  of  the  limit  cycle,  it  can  be  seen  that  any 
counter-example  that  exists  for  the  symmetry  conjecture 
would  have  to  meet  stringent  conditions:   (1)  the  asymmetric 
limit  cycle  (and  its  reflected  cycle)  would  have  to  have 
centroids  on  the  x.  axis,  and  (2)  the  switching  periods 
would  be  unequal.   No  such  counter-example  has  been  found, 
but  no  proof  of  symmetry  exists,  either.   The  question  of 
symmetry  of  friction  limit  cycles  is  therefore  still  open. 

Several  additional  mathematical  tools  can  provide  some 
insight  on  this  problem,  although  none  is  successful  in 
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proving  or  disproving  the  general  case.   Some  notes  about 
these  methods  are  as  follows: 

1.  Degree  Theory:   Since  the  degree  of  an  odd  mapping  is 
odd,  a  possible  attack  on  the  problem  was  to  prove  the 
symmetry  of  limit  cycles  by  proving  uniqueness,  as  follows. 
If  the  limit  cycle  was  not  symmetric,  then  a  reflection 
about  the  origin  would  generate  another  limit  cycle, 
presumably  with  the  same  degree.   Since  the  degree  of  the 
region  containing  the  limit  cycles  must  be  odd  by  odd 
mapping  theorem,  and  the  degree  of  region  is  the  sum  of  the 
degrees  of  the  critical  points  in  the  region,  a 
contradiction  would  be  obtained  since  adding  the  degrees  of 
the  two  limit  cycles  would  give  an  even  answer  (even  +  even 
or  odd  +  odd  are  both  even) .   However,  the  degree  of  any 
limit  cycle  is  actually  zero,  since  the  degree  of  any  point 
on  the  periodic  trajectory  is  zero  (the  Jacobian  of  the 
mapping  must  vanish  since  every  point  on  cycle  is  a  solution 
of  periodic  map;  see  first  part  of  this  chapter) .   The 
oddness  of  the  degree  is  actually  generated  by  the  critical 
point  at  the  origin,  and  limit  cycles  do  not  contribute. 
Therefore,  no  information  about  the  number  or  properties  of 
the  limit  cycles  can  be  obtained  by  considering  the  degree 
of  the  region. 

2.  Fourier  Analysis:   Since  the  even  Fourier  coefficients 
(coefficients  of  even  harmonics)  are  zero  if  and  only  if  the 
function  to  be  expanded  is  half-wave  symmetric,  I  hoped  to 
show  the  symmetry  by  proving  the  even  coefficients  to  be 
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zero  for  the  expansions  of  the  state  trajectories.   This 
approach  became  circular,  however,  since  the  only  way  to 
show  the  even  coefficients  dropped  out  was  to  assume 
half -wave  symmetry. 

Note  that  this  implies  a  possible  weakness  in  the 
describing  function  method,  since  this  method  assumes  at  the 
start  that  the  input  to  the  nonlinearity  is  a  single 
half-wave  symmetric  sinusoid  (the  odd  symmetry  of  the 
nonlinearity  is  also  used  to  show  that  the  cosine  terms  are 
zero  (see  Vidyasagar,  1978)).   An  asymmetric  limit  cycle 
solution  has  unequal  switching  periods,  as  shown  above. 
Therefore,  this  type  of  limit  cycle  might  not  be  predicted 
at  all  by  DF  methods!   An  approach  including  the  second 
harmonic  might  work,  however,  to  better  approximate  the 
asymmetric  input  function. 

3.  Invariant  Surfaces:   By  virtue  of  the  theorem  proved 
above  for  the  2D  case,  if  it  could  be  shown  that  the  limit 
cycle  and  the  origin  (equilibrium  point)  were  both  contained 
in  a  2D  invariant  surface  in  n-dimensional  space,  the 
problem  could  be  reduced  to  the  2D  case  by  a  suitable 
coordinate  transformation  and  solved.   Of  course,  the  limit 
cycle  itself  is  an  invariant  curve,  and  any  closed  cui-ve  is 
topologically  equivalent  to  a  unit  circle  in  the  plane,  but 
it  is  not  known  how  to  construct  such  an  invariant  surface. 
I  examined  the  requirements  for  a  trajectory  to  be  planar 
for  n  =  3,  since  this  is  the  simplest  case  of  an  invariant 
surface,  but  the  case  is  too  specific.   Unfortunately,  a 
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counter-example  already  has  been  presented  for  this  case  in 
n  dimensions  (Example  V-3) ,  although  the  method  may  work  for 
more  specific  cases. 
Summary  of  Results  on  Symmetry 

To  summarize  the  results  of  this  section,  then,  the 
general  case  of  n-dimensional  odd  systems  can  in  fact 
exhibit  asymmetric  limit  cycles.   For  the  more  restrictive 
case  of  odd  Lure  systems  (with  equilibrium  points  at  or  near 
the  origin) ,  no  counter-example  has  yet  been  found  to  the 
conjecture  that  systems  of  this  type  (with  suitable 
conditions  to  reject  trivial  cases)  have  symmetric  limit 
cycles.   Counter-examples,  although  somewhat  contrived, 
might  be  constructed  for  a  system  with  eigenvalues  that  have 
zero  real  parts.   Unfortunately,  the  lack  of  proof  for  the 
conjecture  may  require  the  more  tedious  asymmetric  form  of 
the  piecewise  linear  limit  cycle  equations  be  used — 
otherwise,  a  system  with  an  asymmetric  limit  cycle  might 
incorrectly  be  declared  free  of  limit  cycles  by  using  the 
symmetric  form  of  the  equations. 


CHAPTER  VI 
CONCLUSIONS  AND  RECOMMENDATIONS  FOR  FURTHER  RESEARCH 

Results  and  Conclusions  of  Current  Research 

In  this  section,  a  suiiunary  of  the  results  obtained 
during  this  research  effort  is  provided.   The  summary  also 
identifies,  to  the  best  of  my  knowledge,  those  results 
original  to  this  research. 

It  is  felt  that  these  results  break  new  ground,  since 
the  area  of  study  has  received  little  or  no  attention  by 
researchers  in  control  theory.   While  considerable  progress 
has  been  made  in  opening  up  this  new  area,  and  the  results 
obtained  so  far  are  useful  in  the  practical  design  of 
servomechanisms,  there  exist  significant  open  areas  suitable 
for  further  research. 

The  problem  of  existence  conditions  for  nonlinear 
friction  limit  cycles  was  addressed  first.   In  studying  this 
problem,  the  following  original  developments  were 
accomplished: 

1)  An  eguivalent  control  canonical  model  was  developed 
for  the  general  servomechanism  model  with  friction.   The 
canonical  model  provided  certain  advantages  in  the 
derivation,  and  allowed  the  use  of  a  (well-known)  trick  to 
simplify  the  piecewise  linear  trajectory  equations  (although 
the  control  form  was  not  required  to  make  use  of  this 
simplification)  . 
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2)  Using  this  model,  nonlinear  algebraic  equations 
describing  friction  limit  cycle  trajectories  exactly  were 
constructed  using  the  piecewise  linear  method. 

3)  These  equations  were  shown  to  be  necessary 
conditions  for  the  existence  of  a  (simple)  limit  cycle 
(Theorem  4-1) . 

4)  The  property  of  consistency  was  defined,  which  holds 
when  the  solution  trajectory  is  actually  within  the  proper 
piecewise  linear  region  during  each  switching  period. 

5)  It  was  shown  (Theorem  4-2)  that  a  consistent 
solution  to  the  nonlinear  algebraic  equations  of  Theorem  4-1 
is  a  necessary  and  sufficient  condition  for  the  existence  of 
a  friction  limit  cycle. 

6)  The  method  was  applied  to  a  3D  example  system,  which 
was  shown  to  have  two  friction  limit  cycles,  one  of  each 
standard  type  (pure  sliding  and  sticking) .   The  method  was 
confirmed  by  simulation,  which  exhibited  the  limit  cycle 
predicted  (exactly)  by  the  existence  conditions. 

In  addition  to  the  existence  question,  two  other 
properties  of  the  friction  limit  cycle  were  examined. 
Chapter  V  studied  the  orbital  stability  of  the  limit  cycles 
(i.e.,  whether  trajectories  near  the  orbit  converge  to  it  or 
diverge) ,  and  their  symmetry.   On  the  subject  of  stability, 
the  following  (original)  results  were  obtained: 

7)  An  initial  small  deviation  from  the  limit  cycle  was 
assumed  and  the  resulting  trajectory  for  one  complete  cycle 
was  calculated  using  the  piecewise  linear  technique.   The 
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resultant  deviation  after  one  cycle  was  expressed  as  a 
function  of  the  initial  deviation,  defining  the  so-called 
Poincare  map  (a  standard  method  in  the  literature,  no  known 
application  to  friction) . 

8)  A  linearization  of  the  Poincare  map  was  performed, 
and  the  resulting  matrix  identified  as  the  stability  matrix, 
whose  eigenvalues  determine  the  orbital  stability  of  the 
limit  cycle.   The  calculations  were  performed  for  both  the 
case  of  a  pure  sliding  limit  cycle  and  for  the  sticking 
case;  no  symmetry  property  was  required  or  assumed. 

9)  Although  existing  results  in  a  reference  applied 
only  to  smooth  system  functions,  it  was  proved  that  the 
matrix  derived  by  linearization  of  the  Poincare  map  in  fact 
defines  the  orbital  stability  of  the  friction  limit  cycle. 
An  argument  was  used  based  on  the  complete  map  being  a 
composition  of  smooth  trajectory  pieces,  so  that  convergence 
could  be  proved.   It  seems  likely  that  this  technique  is  a 
standard  approach  to  these  problems,  so  the  author's  proof 
of  orbital  stability  for  the  piecewise  linear  case  is 
probably  not  original,  or  it  is  obvious. 

10)  The  stability  matrix  for  the  3D  example  system  was 
derived,  and  shown  by  simulation  to  accurately  reflect  the 
behavior  of  trajectories  near  the  limit  cycles  in  this 
system. 

The  question  of  the  symmetry  of  limit  cycles  for  odd 
systems  in  general,  and  for  friction  systems,  was  studied. 
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A  number  of  interesting  results  were  obtained  on  this 
question,  which  are  original,  to  the  author's  knowledge: 

11)  Systems  with  odd  system  function  f(.)  have  phase 
portraits  which  are  symmetric  about  the  origin  (a  known 
result) . 

12)  However,  counterexamples  were  identified  which 
showed  that  periodic  orbits  were  not  necessarily  symmetric 
about  the  origin  in  odd  systems. 

13)  Some  conjectures  about  symmetry  in  these  systems 
were  proposed  and  shown  to  be  false:   a)  periodic 
trajectories  were  always  symmetric  about  an  equilibrium 
point,  b)  uniqueness  of  equilibrium  point  implies  symmetry 
about  that  point,  c)  Lure  systems  always  have  symmetric 
limit  cycles,  and  d)  controllable  systems  (in  the  linear 
part)  have  symmetric  orbits.   The  author  is  indebted  to  V. 
M.  Popov  for  an  interesting  counter-example. 

14)  Some  special  cases  were  identified  in  which 
symmetry  held:   a)  if  the  limit  cycle  is  unique,  b)  if  it  is 
unique  in  a  symmetric  region,  or  c)  in  the  2D  case. 

15)  Some  interesting  properties  of  odd  systems  were 
derived:   a)  the  property  of  origin  symmetry  is  equivalent 
to  half-wave  symmetry,  b)  orbits  have  their  centroids  and 
all  equilibrium  points  along  the  Xj^-axis  (for  controllable 
Lure-type  systems) ,  and  c)  symmetry  is  equivalent  to  equal 
switching  periods  (for  friction  cycles). 

The  results  above  were  the  product  of  a  considerable 
amount  of  study  on  this  problem;  unfortunately,  it  was 
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unsuccessful  in  either  proving  symmetry  for  cases  of 
practical  interest  such  as  piecewise  linear  nonlinearities, 
or  in  developing  a  counter-example. 

To  summarize  the  results  of  this  research,  the 
existence  and  certain  properties  of  friction  limit  cycles 
have  been  characterized  by  the  piecewise  linear  method, 
which  provides  exact  results  on  period,  amplitude, 
trajectory,  etc.   The  model  system  is  a  fairly  general 
representation  of  a  servomechanism,  and  can  be  easily 
extended  to  apply  to  related  systems.   Properties  of  the 
limit  cycles  such  as  local  stability  and  symmetry  were 
examined.   Although  some  questions  remain  unanswered  (see 
below) ,  the  results  are  felt  to  be  a  reasonably  complete 
characterization  of  these  limit  cycles,  and  therefore 
provide  an  extremely  useful  tool  in  the  analysis  of  servo 
systems. 

Recommendations  for  Further  Research 

There  are  a  number  of  interesting  questions  in  this 
area  that  are  open.   A  few  of  the  more  interesting  are 
provided  as  suggestions  for  future  researchers: 

1)  The  method  obtained  is  limited  as  a  design  tool 
because  no  closed-form  solution  is  known  for  the  nonlinear 
algebraic  equations.   If  none  can  be  found,  information 
about  solutions  such  as  approximations,  uniqueness,  etc. 
would  be  useful. 

2)  Conditions  on  the  system  that  guarantee  stability 
could  be  found.   For  example,  is  an  unstable  complex  pole 
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pair  required  for  the  existence  of  an  unstable  sliding  limit 
cycle? 

3)  A  proof  that  all  friction  limit  cycles  are  simple, 
or  a  counter-example,  would  answer  the  question  of  whether 
the  existence  conditions  completely  characterize  all  the 
limit  cycles  of  a  system. 

4)  Are  all  sticking  limit  cycles  stable  and  sliding 
ones  unstable? 

5)  What  is  the  relationship  between  the  predictions  of 
other  method  of  nonlinear  analysis  (describing  functions, 
absolute  stability,  etc.)  and  this  method? 

6)  Are  the  necessary  and  sufficient  conditions  stated 
in  Chapter  IV  minimal .  or  are  there  redundancies?   Are  there 
simpler  tests  that  still  form  a  complete  set  of  conditions? 

Two  open  problems  on  the  question  of  symmetry  were 
identified,  requiring  either  a  proof  or  a  counter-example: 

7)  Is  the  stated  conjecture  true,  that  controllable  and 
observable  Lure-type  odd  systems  with  unique  equilibrium 
points  always  have  symmetric  limit  cycles? 

8)  As  a  special  case  of  the  above,  are  friction  limit 
cycles  always  symmetric? 

Indeed,  the  general  question  of  an  exact 
characterization  of  when  odd  systems  have  symmetric  limit 
cycles  remains  open. 

Finally,  the  application  of  this  method  to  friction 
limit  cycles  could  be  duplicated  for  other  piecewise  linear 
nonlinearities  -  a  large  number  exist  that  must  be  accounted 
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for  in  the  process  of  practical  servo  design.   Results  on 
combinations  of  nonlinearities,  which  is  the  typical  case  in 
practice,  would  also  be  highly  useful. 


APPENDIX 
LISTINGS  OF  ANALYSIS  COMPUTER  PROGRAMS 

The  listings  of  two  computer  programs  are  provided  that 
may  be  useful  in  the  study  of  friction  limit  cycles. 
Listing  1  is  a  FORTRAN  program,  using  Newton's  method  of 
solving  nonlinear  equations,  to  calculate  the  switching 
periods  from  the  Chapter  IV  necessary  conditions.   This 
program  has  the  capability  for  solving  the  sliding  case  only 
(no  sticking  assumed) .   A  description  of  the  system  must  be 
provided  in  a  data  file,  and  an  initial  guess  close  to  the 
actual  switching  periods  provided  to  assure  convergence. 

The  second  listing  is  the  FORTRAN  code  describing  the 
state  equations  for  the  3D  example  referred  to  in  the  text. 
Also  included  is  the  subroutine  defining  the  friction  model 
used  in  this  research.   These  subroutines  were  used  with  a 
standard  simulation  package,  including  a  Runge-Kutta  fourth- 
order  integration  routine,  to  simulate  the  example  systems. 
Also  given  is  a  typical  data  file,  showing  model  constants, 
initial  conditions,  and  integration  time  parameters. 
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LISTING  1 
NONLINEAR  EQUATION  SOLVER  FOR  SLIDING  CASE 

C  **************************************************** 

C      PROGRAM  FOR  FINDING  ZEROS  OF  NONLINEAR  FUNCTION  IN 

C  VARS . 

C         (SET  UP  TO  FIND  SWITCHING  PERIODS  OF  SLIDING 

C  FRICTION  LIMIT  CYCLE) 

C 

DIMENSION  XZERO(IO) ,XSTAR(10) ,XGAMM(10) 
COMMON/ 10/ ITERM, IDUMP, IDATA 
COMMON/XDATA/X0(10) ,XT1(10) 

COMMON/CON/CON(10) ,A(100) ,XGAMMA(10) , XL ( 10) , XI ( 100) 
C 

ITERM=12 
IDUMP=14 
IDATA=15 
C 

C  ***  GET  DATA  DEFINING  SYSTEM  FROM  DATA  FILE 
CALL  INPUT (NORDER) 
5     WRITE( ITERM, ' ("  ENTER  START  VALUES  FOR  SWITCHING  ", 
1   "PERIODS") ' ) 
READ(ITERM,)  XZERO ( 1) , XZERO (2 ) 
WRITE (IDUMP, ' ("   START  VALUES  =", 
1    1P2E12.5)')  XZERO(l) ,XZER0(2) 
IERR=0 

CALL  VRNEWT(XZER0,XSTAR,IERR,2) 
IF(IERR.NE.O)  THEN 

WRITE (ITERM, ' ("  lERR  =",I3,";  TRY  NEW  INITIAL", 
1      "  POINT") •)  lERR 
GO  TO  5 
ENDIF 

WRITE (ITERM, ' ("  lERR,  Tl,  T3:")') 
WRITE( IDUMP, ' ("  lERR,  Tl,  T3:")') 
WRITE(ITERM, )  lERR, (XSTAR ( I) , 1=1 , 2) 
WRITE(IDUMP, )  lERR, (XSTAR(I) ,1=1,2) 
C 

C  ***  WRITE  OUT  SWITCHING  STATES 

C         (NOTE  MUST  CONVERT  BACK  TO  PHYSICAL  VARIABLES) 
YN=0. 
YN1=0. 
YNT1=0. 
YN1T1=0. 
DO  9  I=1,N0RDER-1 

XGAMM ( I ) =XGAMMA ( I + 1 ) 
9   CONTINUE 

XGAMM (NORDER) =0. 
DO  10  1=1, NORDER 

YN=YN+XGAMMA ( I ) *X0 ( I ) 
YNT1=YNT1+XGAMMA ( I ) *XT1 ( I ) 
YN1=YN1+XGAMM(I) *XO(I) 
YN1T1=YN1T1+XGAMM ( I ) *XT1 ( I ) 
10   CONTINUE 

WRITE (ITERM, • ("  N,  XO,  XTl:")') 
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99 
98 


WRITE (IDUMP, ' (" 
WRITE (ITERM, 98) 
WRITE (ITERM, 99) 
WRITE (IDUMP, 98) 
WRITE (IDUMP, 99) 
FORMAT (5 (IX, F8. 


N,  XO,  XTl:") ') 

NORDER, (XO(I) ,I=l,N0RDER-2) ,YN1,YN 
(XTl(I) ,I=l,N0RDER-2) ,YN1T1,YNT1 
NORDER, (XO(I) ,I=l,N0RDER-2) ,YN1,YN 
(XTl(I) ,I=l,N0RDER-2) ,YN1T1,YNT1 
4)/3X,5(lX,F8.4) ) 


F0RMAT(I3/5(1X,F8.4)/3X,5(1X,F8.4) ) 
STOP 
END 
*************************************************** 
SUBROUTINE  INPUT (NORDER) 

COMMON/CON/CON(10) ,A(100) ,XGAMMA(10) ,XL(10) , XI (100) 
COMMON/ DIMENS/  NA ( 2 ) , NB ( 2 ) , NC ( 2 ) , ND ( 2 ) , N , NG ( 2 ) 
COMMON/ 10/ ITERM, IDUMP, I DATA 


READ (I DATA,) 
NORDER=N 
READ (I DATA, ) 
READ (I DATA, ) 
READ (I DATA, ) 


C 
C 

c 
c 


N 


(XGAMMA(I) ,I=1,N) 

(XL(I),I=1,N) 

((A(N*(J-1)+I) ,J=1,N) ,I=1,N) 


*** 


ALSO  NEED  TO  SET  UP  DIMENSIONS  FOR  MATRIX  ROUTINES 
NOTE:  ORACLS  ROUTINES  REQUIRE  MATRICES  PACKED  BY 
COLUMNS  IN  ONE-DIMENSIONAL  ARRAY 


NA(1 
NA(2 
NB(1 
NB(2 
NC(1 
NC(2 
ND(1 
ND(2 
NG(1 
NG(2 


=N 
=N 
=N 
=N 
=N 
=N 
=N 
=1 
=1 
=N 


C 

c 

c 

c 


***  INPUT  HOLLERITH  DATA  FOR  TITLE  OF  OUTPUT 
CALL  RDTITL 

***  CREATE  IDENTITY  MATRIX,  ORDER  N 
CALL  UNITY ( XI, NA) 


C* 

c 
c 
c 
c 
c 
c 
c 
c 
c 


RETURN 
END 
********************************************** 

SUBROUTINE  VRNEWT (XZERO, XSTAR, lERR, N) 

NEWTON'S  METHOD  FOR  REAL-VALUED  VECTOR  FUNCTIONS  OF  A 
REAL  VECTOR  VARIABLE 

VARIABLES : 

N      =  DIMENSION  OF  VECTOR  FUNCTION  AND  VARIABLE 

(MAX  =  10) 
XZERO  =  INITIAL  GUESS  AT  SOLUTION 
XSTAR  =  SOLUTION,  CORRECT  TO  6  DECIMAL  PLACES 
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C         lERR   =  0  IF  SOLUTION  FOUND  OK,  ELSE  INDICATES 

C  PROBLEMS : 

C  =  1,2,0R  3  IF  FX  SINGULAR  FOR  SOME  ITERATE 

C  =  4  IF  METHOD  DIVERGES  (DOES  NOT  TERMINATE 

C  IN  MAX  ITERS) 

C  =  11  IF  ORDER  OF  PROBLEM  TOO  LARGE  (N.GT.IO) 

C 

C      SUBROUTINES  NEEDED: 

C         FNDEF  :  DEFINES  FUNCTION  FOR  WHICH  SOLUTION 

C  DESIRED,  F(XSTAR)=0 

C         FXDEF  :  DEFINES  DERIVATIVE  OF  FUNCTION  F  W.R.T.  X 

C         GELIMS:  PERFORMS  GAUSSIAN  ELIMINATION  TO  SOLVE 

C  LINEAR  SYSTEM 

C 

C      STOPPING  CONDITIONS: 

C         MAX  NO.  OF  ITERATION  =  100 

C  DIFFERENCE  BETWEEN  SUCCESSIVE  ITERATES  LESS  THAN 

C  0.5D-7 

C  (NOTE  THAT  THIS  IS  ABSOLUTE  ERROR,  NOT  RELATIVE 

C  ERROR) 

C         F(X)  LESS  THAN  0.5D-11 

C 

DIMENSION  XZERO(IO) ,XSTAR(10) ,XNEW(10) 

DIMENSION  X(10) ,ALPHA(10) 

DIMENSION  BETA(10, 10) ,DELTAX(10) 

COMMON/ 10/ ITERM, IDUMP, IDATA 

DATA  EPSl/0 . 5D-11/ , EPS2/0 . 5D-7/ , ITMAX/100/ 
C 

IF(N.GT. 10) 
IERR=11 
RETURN 

ENDIF 

DO  5  J=1,N 
5     X(J)=XZERO(J) 

1=0 

IERR=0 

WRITE (IDUMP, • ("    I         X(J),J=1,N  ' 
1    "  FJ(X) ,J=1,N    ") ') 
1     CALL  FNDEF (X, ALPHA) 

WRITE (I DUMP,*)  I, (X(J) ,J=1,N) , (ALPHA(J) ,J=1,N) 
IF (VRABS (ALPHA, N) .LT.EPSl) 
DO  2  0  J=1,N 
20       XSTAR(J)=X(J) 
RETURN 
ENDIF 

CALL  FXDEF (X, BETA) 

CALL  GELIMS ( BETA , DELTAX , ALPHA , N , lERR) 
IF(IERR.NE.O) 

RETURN 
ENDIF 

DO  2  5  J=1,N 
25    XNEW(J)  =  X(J)  -  DELTAX (J) 
IF ( VRABS ( DELTAX, N) .LT.EPS2) 
DO  30  J=1,N 
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30       XSTAR(J)  =  XNEW(J) 

CALL  FNDEF(XSTAR, ALPHA) 

WRITE (I DUMP,*)  I, (XSTAR(J) ,J=1,N) , (ALPHA(J) ,J=1,N) 

RETURN 
ENDIF 
1=1+1 
IF(I.GT.ITMAX) 

IERR=4 

RETURN 
ENDIF 

DO  4  0  J=1,N 
40    X(J)  =  XNEW(J) 
GO  TO  1 
END 
C*********************************************** 

SUBROUTINE  FNDEF (X, ALPHA) 

DIMENSION  B(IOO) ,C(100) ,D(10) , 
1   DUMMY(500) ,WK(200) ,IPIVOT(10) 

DIMENSION  ALPHA(IO) ,X(10) 

COMMON/EAT/  EATl (100) ,EAT3 (100) ,EAT(100) ,XINV(100) 

COMMON/CON/CON(10) ,A(100) ,XGAMMA(10) ,XL(10) ,XI(100) 

COMMON/ DIMENS/  NA ( 2 ) , NB ( 2 ) , NC ( 2 ) , ND ( 2 ) , N , NG ( 2 ) 

COMMON/XDATA/X0(10) ,XT1(10) 

COMMON/ 10/ I TERM, I DUMP, I DATA 

DATA  IOP/0/ 
C 

T1=X(1) 

T3=X(2) 
C 

C  ***  FORM  EXP(A*T1),  EXP(A*T3),  EXP(A*T) 
C         (NOTE  MATRICES  ALSO  USED  IN  CALCULATION  OF  JACOBIAN) 

CALL  EXPINT ( A , NA , EATl , NB , C , NC , Tl , lOP , DUMMY ) 

CALL  EXPINT ( A , NA , EAT3 , NB , C , NC , T3 , lOP , DUMMY) 

CALL  MULT(EAT1,NB,EAT3,NB,EAT,NC) 
C 

C  ***  CALCULATE  INV (I-EXP (A*T) ) 
C         (NOTE  MATRIX  ALSO  USED  IN  CALCULATION  OF  JACOBIAN) 

CALL  SUBT(XI,NA,EAT,NA,C,NC) 

CALL  EQUATE ( XI, NA,XINV,NC) 

CALL  GELIM (N , N , C , N , XINV , IPIVOT , 0 , WK , lERR) 

IF(IERR.EQ.l)  THEN 

WRITE(ITERM, ' ("  (I-EAT)  MATRIX  SINGULAR")') 
STOP 

ENDIF 
C 

C  ***  CALCULATE  F1(T1,T3)  AND  F2(T1,T3) 
C         (NOTE  VECTORS  XO  AND  XTl  ARE  SAVED  FOR  CALC. 
C  OF  JACOBIAN) 

CALL  SCALE (EAT3,NB,C,NC, 2.0) 

CALL  SUBT(EAT,NC,C,NC,B,NB) 

CALL  ADD(B,NB,XI,NA,C,NC) 

CALL  MULT(C,NA,XL,ND,C,NC) 

CALL  MULT(XINV,NA,C,NC,D,NB) 

CALL  SCALE(D,NB,X0,NC,-1.0) 
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C         (THIS  IS  XO  VECTOR) 

CALL  MULT(XGAMMA,NG,XO,NC,C,NB) 

ALPHA(l)  =  C(l) 
C 

CALL  SUBT(EAT3,NA,XI,NA,C,NC) 

CALL  MULT(C,NC,XL,ND,B,NB) 

CALL  MULT(XINV,NA,B,NB,C,NC) 

CALL  MULT(EAT1,NA,C,NC,D,NB) 

CALL  SCALE (D,NB,C,NC, 2.0) 

CALL  ADD(C,NC,XL,ND,XT1,NB) 
C        (THIS  IS  XTl  VECTOR) 

CALL  MULT (XGAMMA,NG, XTl, NB,C,NC) 

ALPHA (2)  =  C(l) 
C 

RETURN 

END 
C****** ************************ ************* 

SUBROUTINE  FXDEF (X, BETA) 

DIMENSION  X(10) ,BETA(10,10) 

DIMENSION  B(IOO) ,C(100) ,D(10) ,TEMP(10) 

COMMON/EAT/  EATl(lOO) ,EAT3 (100) ,EAT(100) ,XINV(100) 

COMMON/CON/CON (10) ,A(100) ,XGAMMA(10) ,XL(10) ,XI(100) 

COMMON/XDATA/X0(10) , XTl (10) 

COMMON/ DIMENS/  NA (2 ) , NB (2 ) , NC (2 ) , ND (2 ) , N, NG (2 ) 
C 

C  ***  CALCULATE  2X2  JACOBIAN  MATRIX  OF  PARTIAL  DERIVATIVES 
C        OF  Fl  AND  F2  WITH  RESPECT  TO  Tl  AND  T3 
C 

C  ***  D(F1)/D(T1) 
C         (TEMP  HOLDS  X(T1)  -  L) 

CALL  SUBT(XT1,ND,XL,ND,TEMP,NB) 

CALL  MULT (EAT3,NA, TEMP, NB,C,NC) 

CALL  MULT(A,NA,C,NC,B,NB) 

CALL  MULT(XINV,NA,B,NB,D,NC) 

CALL  MULT(XGAMMA,NG,D,NC,C,NB) 

BETA(1,1)  =  C(l) 
C 

C  ***  D(F2)/D(T1) 
C        (D  HOLDS  XINV*A*EAT3*(XT1-L)  ) 

CALL  MULT(EAT1,NA,D,NC,B,NB) 

CALL  MULT(XGAMMA,NG,B,NB,C,NC) 

BETA(2,1)  =  C(l) 

CALL  MULT(A,NA,TEMP,ND,B,NB) 

CALL  MULT(XGAMMA,NG,B,NB,C,NC) 

BETA(2,1)  =  BETA(2,1)  +  C(l) 
C 
C  ***  D(F1)/D(T3) 

CALL  ADD(XT1,ND,XL,ND,B,NB) 

CALL  MULT(EAT3,NA,B,NB,C,NC) 

CALL  MULT(A,NA,C,NC,B,NB) 

CALL  MULT(XINV,NA,B,NB,D,NC) 

CALL  MULT(XGAMMA,NG,D,NC,C,NB) 

BETA (1,2)  =  C(l) 
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C  ***  D(F2)/D(T3) 

C        (D  HOLDS  XINV*A*EAT3*(XT1+L)  ) 

CALL  MULT(EAT1,NA,D,NC,B,NB) 

CALL  MULT(XGAMMA,NG,B,NB,C,NC) 

BETA (2, 2)  =  C(l) 
C 

C  ***  RETURN  MATRIX 
C 

RETURN 

END 
C  ****************************************************** 

SUBROUTINE  GELIMS ( Al , XI , Bl , N , lERR) 

DIMENSION  AORIG(10,10) ,AINV(10,10) ,X(10) , 
1   B(10) ,A(10,10) 

DIMENSION  Al(10,10) ,B1(10) ,X1(10) 

COMMON/MATRIX/ AORIG , A , AINV , B , X 
C 

C   GAUSSIAN  ELIMINATION  FOR  SYSTEMS  OF  ORDER  . LE .  10 
C      Al  =  10X10  DIMENSION  SYSTEM  MATRIX 
C      Bl  =  10X1  DIMENSION  DATA  VECTOR 
C      XI  =  10X1  DIMENSION  SOLUTION  VECTOR 
C      N   =  ACTUAL  ORDER  OF  SYSTEM 

C      lERR  =  ERROR  CODE,  =0  IF  OK,  =11  IF  ORDER  TOO  LARGE, 
C  =1,2,3  IF  DIFFICULTIES  IN  GAUSSIAN 

C  ELIMINATION  (SINGULAR  MATRIX,  ETC.) 

C 

IF(N.GT.IO) 
IERR=11 
GOTO  900 
ENDIF 

DO  10  1=1, N 
B(I)=B1(I) 
DO  10  J=1,N 

A0RIG(I,J)=A1(I,J) 
10    CONTINUE 

CALL  FORWARD (N, I ERR) 
IF(IERR.NE. 0)  GOTO  900 
CALL  BACK(N,IERR) 
IF(IERR.NE.O)  GOTO  900 
DO  20  1=1, N 
X1(I)=X(I) 
2  0    CONTINUE 

RETURN 
900   RETURN 
END 
C  ************************************************ 

SUBROUTINE  FORWARD (ND, lERR) 

DIMENSION  AORIG (10, 10) , AINV (10, 10) ,X(10) , 
1   B(10) ,A(10,10) 
COMMON/MATRIX/AORIG , A, AINV, B, X 
DIMENSION  S(10) 
INTEGER* 3  R 
COMMON/INDEX/R(10) 
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C      INITIALIZE  INDEX  VECTOR  AND  COPY  ORIGINAL  MATRIX 

C  (SO  ELIMINATION  IS  NON-DESTRUCTIVE) 

C 

DO  10  K=1,ND 

DO  11  J=1,ND 
11       A(K,J)=AORIG(K,J) 
10    R(K)=K 
C 

C      FORWARD  ELIMINATION 
C 

DO  100  K=1,ND-1 
C 

C         CALCULATE  ROW  SCALE  FACTORS 
C 

DO  60  I=K,ND 
S(I)=0.0 
DO  55  J=K,ND 

IF(ABS(A(R(I) ,J) ) .GT.S(I) ) 
1  S(I)=ABS(A(R(I) ,J) ) 

55  CONTINUE 

C 

C  CHECK  FOR  SINGULAR  A   (ANY  S(I)  =  0) 

C 

IF(S(I) .LE.l.OE-100) 
IERR=1 
GO  TO  900 
ENDIF 
60       CONTINUE 
C 

C         FIND  PIVOT  ROW 
C 

ISTAR=K 

TPMAX=0 . 

DO  70  I=K,ND 

TPIVOT=ABS ( A (R ( I ) , K) ) /S ( I ) 
IF (TPIVOT . GT . TPMAX) 
TPMAX=TPIVOT 
ISTAR=I 
ENDIF 
7  0        CONTINUE 
C 

C         CHECK  FOR  SINGULAR  A   (TPMAX  =0.) 
C 

IF (TPMAX. LE.l.OE-100) 
IERR=2 
GO  TO  900 
ENDIF 
C 

C         IMPLICIT  PIVOTING  BY  EXCHANGING  INDICES  IN 
C  INDEX  VECTOR 

C 

TEMPER (K) 

R(K)=R(ISTAR) 

R(ISTAR)=TEMP 
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C 

C         CALCULATE  PIVOTS  AND  PERFORM  ROW  ELIMINATION 

C  (KTH  STEP) 

C 

XPIVOT=A(R(K) ,K) 
DO  80  I=K+1,ND 

A(R(I),K)  =  A(R(I) ,K)/XPIVOT 
DO  75  J=K+1,ND 

A(R(I),J)  =A(R(I),J)  -  A(R(I)  ,K)*A(R(K)  ,J) 
75  CONTINUE 

80       CONTINUE 
100   CONTINUE 
C 

C      NORMAL  RETURN 
C 

IERR=0. 
RETURN 
C 

C      ERROR  RETURN 
C 

9  00   RETURN 

END 
C  ************************************** 

SUBROUTINE  BACK (ND, lERR) 

DIMENSION  AORIG(10,10) ,AINV(10,10) ,X(10) , 
1   B(10) ,A(10,10) 
COMMON/MATRIX/ AORIG , A , AINV , B , X 
INTEGER*3  R 
COMMON/ INDEX/ R (10) 
C 

C      TRANSFORM  THE  DATA  VECTOR  USING  PIVOTS  STORED  IN  A 
C 

DO  20  K=1,ND-1 

DO  10  I=K+1,ND 

B(R(I))=B(R(I))  -  A(R(I),K)*B(R(K)) 

10  CONTINUE 
20    CONTINUE 

C 

C      BACKWARDS  SUBSTITUTION  FOR  THE  SOLUTION  X 

C 

IF(ABS(A(R(ND) ,ND) ) .LE.l.OE-100)  GO  TO  900 
X(ND)  =  B(R(ND) )/A(R(ND) ,ND) 
DO  50  K=2,ND 
L=ND+1-K 
X(L)=B(R(L)) 
DO  45  J=L+1,ND 

X(L)=X(L)-A(R(L) ,J)*X(J) 
45       CONTINUE 

IF(ABS(A(R(L) ,L) ) .LE.l.OE-100)  GO  TO  900 
X(L)  =  X(L)/A(R(L),L) 
50    CONTINUE 
IERR=0 
RETURN 
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C      ERROR  EXIT 
C 
900   IERR=3 

RETURN 

END 
C  ***************************************** 

REAL  FUNCTION  VRABS(X,N) 
DIMENSION  X(l) 
VRABS  =0.0 
DO  10  1=1, N 
10    IF(ABS(X(I) ) .GT.VRABS)  VRABS  =  ABS(X(I)) 
RETURN 
END 
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LISTING  2 
EXAMPLE  SIMULATION 


I.  STATE  EQUATIONS; 


SUBROUTINE  EQNS 
C  ******************************************* 

C  *     MODEL  FOR  3D  SYSTEM  WITH  NONLINEAR  FRICTION  FOR 

C  *        LIMIT  CYCLE  STUDIES 

C  ****************************************** 

COMMON  MPTS,NPTS,V(2000) ,ITITLE(24) , KPRNT, KSAV, KPLOT, 

1  NCON,CON(300) , NSTATE, XIN (2 ) ,XERR(2) ,X0UT(2) , 

2  XSPEC(5) ,X(75) ,XD(75) , T, DT, TMAX, DTOUT, IFIRST, IRK 
COMMON/ DIGFIL/TAU , NDSTATE , XDIG (75), YDIG (75) 
COMMON/PGMCON/C(100) ,ISW(80) ,H(2,75) 
COMMON/IOOI/LDI , LDO, LDM, LDMO, LDP 

DATA  XDL,TL/0. ,-1./ 
DATA  IPRT/0/ 
C  ******************************************** 

C  * 

C  *   THIS  SUBROUTINE  CONTAINS  THE  USER-DEFINED  SYSTEM 

C  *        STATE  MODEL 

C  * 

c  *************************************** 

c 

C  ****  SYSTEM  INPUT  DEFINED  HERE 
C 

XIN(1)=0.0 

XIN(2)=0.0 

c 

C  ****  SYSTEM  STATE  MODEL 
C 

XD(1)  =  X(2) 
XD(2)  =  X(3) 

TORQUE  =  -C0N(1)*X(1)  -  CON(2)*X(2)  -  CON(3)*X(3) 
FRICT  =  UNLUBED(1.  ,C0N(4)  ,C0N(5)  ,1.E-3,DT,XDL,X(3)  , 
1   TORQUE, 1,T,TL,ISBIT) 
XD(3)  =  TORQUE  -  FRICT 
C 

C      WRITE(3,*)  (XD(I) ,1=1,3) 
IF(ISW(1) .EQ.l) 
DO  31  J=l,3 

XD(J)  =  -XD(J) 
31    CONTINUE 
ENDIF 
IF(ISW(4) .EQ.l)  THEN 

IF(ABS(X(3) ) .LE.1.E-3.AND.X(2) . GT . 0 . 0 
1  .AND.IPRT.EQ.O)  THEN 

IPRT=1 

WRITE(20,)  T,X(1),X(2) 
ELSEIF(ABS(X(3) ) .GT.l.E-3)  THEN 

IPRT=0 
ENDIF 
ENDIF 
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C 

c  **** 


c  **** 


1000 


OUTPUT, ERROR, OTHER  STANDARD  VARIABLES  DEFINED  HERE 
XERR(1)=  -X(l) 
XERR(2)=  -X(2) 
DO  5  1=1,5 
XSPEC(I)=0.0 
XSPEC(1)=  X(3) 
XSPEC(2)  =  X(3) 
XSPEC(3)  =  FRICT 

OUTPUT  DEFINED  AS  LINEAR  SUM  OF  STATES 
XOUT(1)=0.0 
XOUT(2)=0.0 
DO  1000  I=1,NSTATE 
XOUT(2)=XOUT(2)+H(2,I)*X(I) 
X0UT(1)=X0UT(1)+H(1,I)*X(I) 
RETURN 
END 


C  **** 


II.  FRICTION  MODEL: 


C  ** 


c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


*** 

FUNCTION  UNLUBED  (  J,  CF,  SF,  EPS,  DT,  XDLAST,  XD, 
1  TORQUE,  FLAGO,  TIME,  TMLAST,  STKNBIT  ) 

INTEGER  FLAGO,  STKNBIT 
REAL  J 


USER  NOTES: 

J         =  INERTIA  OF  LOAD/MOTOR  CAUSING  FRICTION. 

CF        =  COULOMB  FRICTION  LEVEL. 

SF        =  STICTION  LEVEL. 

EPS       =  THE  STICTION/ COULOMB  BREAKPOINT. 

DT        =  SIMULATION  INCREMENTAL  TIME  STEP. 

XDLAST    =  USUALLY  INITIALIZED  TO  INITIAL  VALUE 

OF  XD.   CON(X)  IS  A  CONVENIENT  WAY  TO 

INITIALZE,  OR  THE  "CONST"  SUBROUTINE 

COULD  BE  USED  TO  SET  A  VARIABLE  OF  YOUR 

CHOICE  EQUAL  TO  THE  INITIAL  VALUE  OF 

XD. 
XD        =  VELOCITY  OF  LOAD/MOTOR  CAUSING  FRICTION. 
TORQUE    =  INPUT  TORQUE  WHICH  WILL  ACCELERATE 

LOAD/MOTOR  BEFORE  COULOMB  FRICTION 

IS  SUBTRACTED  OUT. 
FLAGO     =  FLAG  BY  WHICH  THE  USER  IMPLEMENTS  THE 

STICTION  FUNCTION. 
FLAGO  =  1  WILL  CAUSE  THE  VELOCITY,  XD, 

TO  BE  SET  EQUAL 

TO  0.0  WHENEVER  THE  VELOCITY  CHANGES 

SIGN. 
STKNBIT   =  STICTION  BIT.   THIS  FLAG  SET  TO  1  WHEN 

IN  THE  STICTION  MODE.   THE  PURPOSE  OF 


135 


C  THIS  FLAG  IS  TO  INSURE  THAT  THE 

C  INTEGRATION  ROUTINE  DOES  NOT  CHANGE  XD 

C  TO  BE  SLIGHTLY 

C  DIFFERENT  FROM  0.0  UNTIL  THE  INPUT 

C  TORQUE  EXCEEDS  THE  STICTION  LEVEL,  SF. 

C 

C 

C 


C 

c 
c 


c 
c 

c 

c 


IF(  XD  .NE.  0.0  ) 

IF(  STKNBIT  .NE.  1  ) 

IF(  XDLAST*XD  .GE.  0.0    .OR.    FLAGO  .EQ.  0  ) 


IF(  ABS(XD)  .GE.  EPS  )  THEN 
UNLUBED  =  SIGN(  CF,  XD  ) 
ELSE 

UNLUBED  =  SIGN(  ABS (XD) * (CF-SF) /EPS  +  SF, 
1  XD  ) 

END  IF 
C 

C  CHECK  IF  ACCELERATION  TORQUE  IS  ZERO  (OR 

C  INSIGNIFICANT) . 

IF(  ABS(TORQUE-UNLUBED)  .LT.  l.OE-12) 
1  GO  TO  10 

C  FIND  TIME  TO  DECELERATE  XD  TO  ZERO. 

TIMETOZ  =  -J*XD/(TORQUE-UNLUBED) 
C  SET  XD  TO  ZERO  IF  TIMETOZ  <  DT/2. 

IF(  TIMETOZ  .LT.  0.0  .OR. 

1  TIMETOZ  .GE.  DT/2.0  .OR.  FLAGO  .EQ.  0  ) 

2  GO  TO  10 


END  IF 

END  IF 

XD  =  0.0 
END  IF 
C 

STKNBIT  =  1 

IF(  ABS (TORQUE)  .GT.  SF  )  THEN 
UNLUBED  =  SIGN(  SF,  TORQUE  ) 
STKNBIT  =  0 
ELSE 

UNLUBED  =  TORQUE 
END  IF 
C 

10  IF(  TIME  .NE.  TMLAST  )  TMLAST  =  TIME 

XDLAST  =  XD 
END  IF 
RETURN 
END 
C  ***** 


III.  TYPICAL  DATA  FILE; 
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FRICTION 

3,0 

1 

0.0,0.001 

100,1000, 

0,1,0.0 

8 

1 

8 

1 
2 

3 

El 

E2 

SI 

S2 

II 

1.0 

0.0 

80 

0  0  0 

0  0 

0  0 

0  0 
1.0 
0.0 
8 


0,100 
200 


,100. 


0 
0 
0 


LIMIT  CYCLES  -  3D,  INTEG.  ; TITLE  (24A3) 

;N0.  OF  STATES,  DISCRETE  STATES  (F.F.) 
; INTEGRATION  ROUTINE  (F.F.;0=PC) 
; TO , DT , TMAX , TAU  (F.F.) 
;KPRNT,KSAV,KPLOT  (F.F.) 
STATISTICS  OPTION,  KSTATS ,  STATS.  DELAY 
NO.  OF  VARIABLES  TO  BE  PRINTED  (F.F.) 
LINEPTR  PLOT  OPTION  ( 1=N0  PLOT,  0=PLOT) 
NO.  OF  VARIABLES  TO  BE  PLOTTED  (F.F.) 


VARIABLES  TO  BE  PRINTED,  PLOTTED 


1. 
1. 


0.20 
0.0 

0    0 

0    0 

0    0 

0    0 
0.0 
0.0 

0 
2 


0  0 

0  0 

0  0 

0  0 


0.0  0.0 

0.0  0.0 

;N0.    OF    SWITCHES     (F, 
0000000000 
0    0    0    0    0 


0  0  0 
0  0  0 
0  0  0 
0.0 
0.0 
;  NO. 
0 
0 


0 
0 
0 


0  0  0  0  0 
0  0  0  0  0 
0.0 
0.0 
OF  CONSTANTS  (F 
0.9        1.0 
0.0        0.0 


) 


INIT  CONDS  4E10.3 
ICS  -  DISC  STATES 

SWITCH  DEFS  (2012) 


OUTPUT  VECTOR  1 
OUTPUT  VECTOR   2 

F.) 

;  MODEL  CONSTANTS 
;  8 
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